Adhesion governs to a large extent the mechanical interaction between a cell and its microenvironment. As initial cell spreading is purely adhesion driven, understanding this phenomenon leads to profound insight in both cell adhesion and cell-substrate interaction. It has been found that across a wide variety of cell types, initial spreading behavior universally follows the same power laws. The simplest cell type providing this scaling of the radius of the spreading area with time are modified red blood cells (RBCs), whose elastic responses are well characterized. Using a mechanistic description of the contact interaction between a cell and its substrate in combination with a deformable RBC model, we are now able to investigate in detail the mechanisms behind this universal power law. The presented model suggests that the initial slope of the spreading curve with time results from a purely geometrical effect facilitated mainly by dissipation upon contact. Later on, the spreading rate decreases due to increasing tension and dissipation in the cell' s cortex as the cell spreads more and more. To reproduce this observed initial spreading, no irreversible deformations are required. Since the model created in this effort is extensible to more complex cell types and can cope with arbitrarily shaped, smooth mechanical microenvironments of the cells, it can be useful for a wide range of investigations where forces at the cell boundary play a decisive role.
The dynamics of initial cell spreading – that is during the first few minutes – are governed by energy release through binding events of cell surface molecules, rather than by active cellular processes such as e. g. tension generated by stress fibers. These molecular binding events dominate the total adhesion energy of the cell. This adhesion creates a pulling effect that in turn generates strong local forces which result in deformations of the actin cortex. The dynamics of initial cell spreading (the increase of radius of the contact area with time) universally correspond to an early (), and a later () power law behavior [1]. It is only at an advanced stage when the cell is already moderately spread out that active pulling of actin stress fibers on focal adhesion complexes will reinforce cell spreading, depending on the cell type in question, see e. g. [2]. The viscoelastic behavior of the cell boundary is determined not so much by the cell membrane itself but by the intracellular cytoskeleton, or, in the case of red blood cells (RBCs), a network of spectrin filaments directly underlying the membrane [3], [4]. A model that can be used for describing cellular mechanics should be able to accurately describe the mechanical interactions that take place at the cell boundary, i. e. the contact interface with its substrate, the extracellular matrix or surrounding cells. Lattice-free, particle-based methods can describe the interaction forces and the resulting movement and deformation of particles in a natural way. At a point of contact between two particles, contact forces are calculated explicitly based on an appropriate contact force model. From these forces, movement of the particles is calculated by integrating the equation of motion. In the simplest approach, particles are assumed to be spherical. In that case, contact forces can be directly calculated from the sphere-sphere overlap distance (are the radii of the spheres and the spacial coordinates of their centers). Calculating contact forces for non-spherical shapes is more challenging: approximations have to be made for the contact force model and it is not trivial to calculate a meaningful overlap distance for all cases. Arbitrary shapes have been modeled by using combinations of connected overlapping spheres [5] or by using polyhedra or poly-arcs, and calculating a contact force proportional to the overlapping volume of the shapes [6], [7]. Besides, the surface of an arbitrary shape can be approximated by sampling points [8]. For each sampling point, a contact force can be calculated based on the indentation in the surface of another object. Disadvantages of using sampling points include that it is hard to directly compare it to a physical contact model such as the Hertz model for spheres, that they generally do not allow to reach complete force equilibrium, and that the precision of the approximation of the contact depends crucially on the local density of nodes, so that the contact parameters need to be re-scaled for different node densities [8]. We present a novel computational framework for describing the mechanical behavior of cells with an emphasis on the interaction between the cell and its environment. Although we only apply this model to cell spreading on a flat surface, the current implementation already allows for more complex settings of interaction with arbitrarily shaped smooth bodies, and cell-cell interaction. The main novelty of the method developed in this work lies in the fact that we calculate contact between a triangulated surface with “rounded” triangles reflecting the local curvature of a cell and its microenvironment by applying Maugis-Dugdale theory (see section “Maugis-Dugdale theory”) to all contacting triangles. To apply this adhesive contact model for the triangulated surfaces in our models, we build on the following six ideas (see section “Contact mechanics of a triangulated surface”): This novel contact model is combined with a new implementation derived from existing mechanical models for red blood cells, mainly from Fedosov et al. [3], [10]. That model has been previously computed using a dissipative particle dynamics (DPD) solver, a different meso-scale simulation method. The mechanical model of the cortex of the RBC includes finitely extensible nonlinear elastic (FENE) connections and viscous dissipation between the nodes of the triangulation, volume conservation and surface area conservation, as well as bending resistance – see section “Elastic model of the cortex”. Finally, we apply this newly developed method to an in-depth computational investigation of RBC spreading (see Figure 1 and supplementary Videos S1 and S2) as reported by both Hategan et al. [11] and Cuvelier et al. [1] in order to unravel the governing mechanisms.
In the experiments reported by Cuvelier et al. [1], biotinylated RBCs were osmotically swollen to become spherical and the change of the radius of the contact area with time was measured for spreading on a streptavidin coated surface. To compare to the spreading dynamics reported in that paper as well as by Hategan et al. [11] (where the cells spread on a polylysine coated surface), we set up simulations of the described model with the parameters as given in table 1. The red blood cell is modeled with a viscoelastic cortex including bending stiffness and Maugis-Dugdale contact interactions. Most parameters in table 1 are taken directly from the literature as indicated. The effective range of interaction (see equation 8) was estimated at by interpolating from [14] for cells with a radius of ≈3 µm. The cortex Young' s modulus used in the Maugis-Dugdale model is the material stiffness of the phospholipid-spectrin complex (the elasticity of the deforming membrane is already taken into account by the FENE potentials). This material stiffness can be assumed to be much higher compared to the whole cell' s Young' s modulus and is set at a value of. The parameters for the cortex are validated by performing the cell stretching and relaxation experiments explained in the previous section “Validation of the RBC cortex model”. A view on three stages of the cell spreading for both biconcave and sphered RBCs is presented in Figure 1. Note that the volume of the biconcave RBC is only about of the volume of the sphered RBC. As a result of that, for the sphered RBC, the final height of the spread-out cell is greater and it has a higher angle of contact compared to the final shape of the initially biconcave RBC. For this simulation, a triangulation based on a five-fold subdivision of an icosahedron was used – see section “Generating triangulated meshes of cells”. This level of mesh refinement is required to reproduce the final high curvatures at the edge of the contact area when the cell is fully spread out: The triangles at the edge have encompassing spheres with radii of ca. , while Hategan et al. [15] report a typical radius of the rim for this situation of, which is of comparable order of magnitude. The shape of the final spread-out cell is a spherical cap. By fitting a sphere through the top % of the nodes, the effective contact angle [16] can be estimated. For the modeled RBC, we calculate an effective contact angle of, which corresponds reasonably well to the measured effective contact angle of around 60° [15]. Figure 4 shows the power-law behavior of the sphered RBC spreading in double logarithmic representation. The “contact radius” of the RBC in these and the following figures is calculated from the sum of all the triangles' areas which are in contact by defining. The spreading dynamics of the model match the experimentally observed cell spreading [11] very well. Figure 5 summarizes the influence of varying one parameter at a time for the most influential parameters of the model starting from the base parameter set reported in table 1. Its first sub-figure (a) shows simulation results of cell spreading for different values of the cell-substrate adhesion strength. A lower adhesion strength results in a lower final contact radius, but also makes the spreading slower. However, the power law behavior as reported by Cuvelier et al. [1] stays well conserved for different adhesion strengths. The influence of the FENE stretching constant is shown in Figure 5 (b). In the range of the RBC FENE constant (in the order of 1 µN/m), the influence of on the spreading dynamics is comparatively small. For larger deviations, higher values of limit the final spreading radius to a lower value, or conversely, lower values allow the cell to spread considerably more. A FENE connection is also characterized by the maximal stretch (Figure 5 (c) ), which expresses the maximal extension of the spring, at which the FENE force diverges (equation 22). The initial spreading dynamics are not affected by the precise value of, but the final spreading radius is. For higher values of, the same tension in the cortex corresponds to a larger extension and therefore a larger radius of the spread out cell. The effect of the bending stiffness on RBC spreading is shown in Figure 5 (d). A higher bending resistance of the cortex speeds up cell spreading, the probable reason being that, through resisting to bending, the cortex keeps the contact angle within the effective range of adhesive interaction close to 180°. This range is of the order of for microscopic biomolecular surfaces [14]. It should be noted that for a theoretical vesicle with bending resistance, the actual contact angle is always 180° [16]. However, for a real RBC, the width of the adhesive spreading front is non-zero and determined by the effective range of interaction. This effective adhesive range is taken into account in Maugis-Dugdale theory (equation 8) and relates the maximal adhesive tension at the edge of contact to the total work of adhesion. The normal friction coefficient is determined by the energy dissipation when adhesive contact is initiated. The dissipation is caused by snap-in-contact events when adhesion molecules form bonds, and the hysteresis arising from unbinding stochastically again [14]. In Figure 5 (e), the effect of changing on the RBC spreading dynamics is shown. As could be expected, a lower value of diminishes the energy dissipation due to adhesion and therefore increases the rate of cell spreading. However it does not change the initial power law behavior of cell spreading. Finally, in Figure 5 (f), the effect of the local area constraint on the spreading dynamics is shown. When the value of is too low, degenerate triangle shapes can arise with a strongly decreased area. This will result in an underestimation of the final spreading radius. It can be observed that for values of, the local area of the triangles is sufficiently well conserved and the predicted spreading dynamics are not affected. In Figure 6 (a), the outward normal pressure on the nodes is visualized for three distinct phases of the cell spreading process for a sphered RBC. The normal pressure is defined here as the magnitude of the sum of all conservative forces (on the left-hand side in the equation of motion, 31) in the nodes projected onto the normal in that node – therefore this normal pressure is dominated by contact forces, where adhesive ones yield a positive (outward) pressure in this case. Figure 6 (b) shows the in-plane tension (in) of the cortex (further denoted as cortex tension, and not to be confused by the adhesive tension, given by Maugis-Dugdale theory, see equation 7) at the same time points. This tension is characterized by the FENE force at the inter nodal connections. Positive forces in these connections correspond to tensile stress in the cortex, while negative values are associated with compressive stress: (1) where is the number of FENE connections of node and is the inter-nodal distance (see e. g. [17]). At the spreading dynamics correspond to the power law regime. At this stage, adhesive forces are strong especially at the edge of contact, but also in the entire rapidly increasing circular contact area. The elastic energy stored in the membrane at this point in time is very low, as the stretch and bending in the membrane is small. As a result, almost all the energy dissipation (see section”Equation of motion”) takes place in the contact area. At, a distinct adhesive edge can be observed, in which the magnitude of forces is much stronger than in the inner circle of the contact area, where the contact potential is already nearly minimal. At the edge, the cortex' s bending stiffness provides resistance to the strong adhesive tension. Meanwhile, the upper spherical cap is being stretched while at the plane of contact the membrane – together with the substrate it is adhering to – is under compressive stress. At this stage, energy dissipation takes place not only at the substrate interface, but also in the entire stressed cortex. As a result of this, the spreading slows down to a lower rate than the power law regime. At, spreading has stopped and the cell has reached equilibrium. The forces at the nodes are zero, and the adhesive tension at the edge of contact is being balanced out by the elastic stress in the RBC membrane/cortex. The cortex in the spherical cap is under strong tensile stress and the stretch in the connections is close to its maximal value. At the substrate interface, compressive stresses have built up even more. For an elastic substrate, these compressive forces will cause radial inwards deformation of the substrate, as has been observed in traction force microscopy measurements [18], [19] – although these experiments concern late cell spreading. It should be noted that the maximal normal pressure at the nodes – occurring in the first stage of cell spreading – corresponds to a force in the order of, which is in the range of the force applied in the stretching simulations which were used to validate the model parameters of the elastic cortex, see section “RBC stretching experiments”.
First, with regard to the performance of the newly developed model for a triangulated, deformable cell obeying Maugis-Dugdale contact tractions, we conclude that: The modeling technique described in this work has a number of limitations: Finally, regarding the initial dynamics of cell spreading, we find: Although the model as shown is restricted to RBC spreading dynamics, we expect that these conclusions can be generalized to other cell types: the same key mechanical components are present in other systems as well, and despite the fact that other cells' cytoskeletons are more complex and the cells can dissipate energy through “active biological processes”, we expect the initial cell spreading phase to be still characterized by contact dissipation. Eventually, stress in the membrane/cortex will build up as well and through this, the cell will dissipate energy in the entire cortical shell. However, it is possible that this dissipation involves irreversible deformation in the cortex.
For two spherical asperities in contact or one asperity in contact with a flat surface (see Figure 7), Maugis-Dugdale (MD) theory can be used to describe the contact mechanics [24]. This theory captures the full range between the Derjaguin-Muller-Toporov (DMT) zone of long reaching adhesive forces and small adhesive deformations to the Johnson-Kendall-Roberts (JKR) limit of short interaction ranges and comparatively large adhesive deformations in the transition parameter. This transition parameter relates to the Tabor coefficient by a factor of [25]. (3) In equation 3, is the maximum adhesive tension (measured in Pa) from a Lennard-Jones potential, (in J/m2) the adhesion energy, is the reduced radius of the asperities and the combined elastic modulus: (4) The (repulsive) Hertz pressure associated with a contact of radius (see Figure 7) is given by (5) Assuming a spherical asperity – and therefore a circular contact area – the total Hertz force can be calculated by integrating equation 5 over the complete circular contact area with radius, which yields the total Hertz force: (6) An adhesive stress can be formulated as [24], [26]: (7) In the Maugis-Dugdale model, local adhesion tension is assumed to be independent of the overlap until a cut-off distance. If the asperity is further than away from the flat surface, the adhesive tension drops to zero. Therefore, is related to the adhesion energy as: (8) is the total work of adhesion, i. e. the work required to move the asperity away from the surface and out of contact. To pull a small area out of contact, the required work is: (9) The total (global) adhesive force is the integral over the adhesive zone with radius (see Figure 7), which according to [26] becomes: (10) The force in equation 10 is dependent on. As equation 10 expresses the global adhesive force of the complete asperity, it is not a constant force, but through dependent on the indentation. To calculate the adhesive radius from the actual geometrical contact area with radius, the height at the edge of the adhesive zone can be used. Substituting both repulsive and adhesive pressures at (see equation 5 and 7) this yields [25]: (11) where. In general, to calculate both and from a given state of the contact, one needs to solve this equation simultaneously with the equation for the net contact force [25]: (12) A very well validated contact model for soft, adhesive bodies like cells, the JKR theory [27]–[29], is a limiting case of Maugis-Dugdale theory for negligible cutoff-distance for the adhesive interaction (or). It has therefore a parameter less than MD theory. The adhesive pressure according to JKR (compare to equation 7) is (13) Note that this pressure diverges at. Summarizing the Maugis-Dugdale theory for an adhesive contact, one considers three distinct zones: The meshes used in this work are derived from spherical shapes by subdividing an icosahedron and projecting the nodes on a sphere [30]. In a subdivision, each triangle gets split into four triangles as is illustrated in Figure 8. Here it is shown how one triangle with an encompassing sphere matching the local curvature of the cell, is split into four triangles. Since the local curvature is kept, the new triangle nodes are all located on the surface of the same encompassing sphere. Every subdivision of an icosahedron has only twelve nodes with a five-fold connectivity and slightly longer distances to their neighbors; otherwise, the mesh is perfectly regular with six-fold connectivity and is ideal for curvature calculations (see section “Local curvature of the 3D shape”) as reported by [31]. The bi-concave shape of an RBC can be obtained by reducing the volume of the sphere to approximately 60%, and letting a system of linear springs with appropriately chosen parameters relax again. This is the reverse process to the well described technique of RBC sphering, see e. g. [32]. We use meshes of either four or five subdivisions of an icosahedron, corresponding to 642 and 2562 nodes, respectively. In the deformable cell model, the cortex nodes interact through viscoelastic potentials. In the most simple approach, a linear elastic spring could be used. For a given displacement of nodes and, the elastic spring force over a connection is: (20) in which and are the actual distance and equilibrium distance between connected node and. The linear spring stiffness is called. For red blood cells, two non-linear spring models have been used in literature: the finite extensible non-linear elastic model (FENE) and the worm-like chain model (WLC) [3]. These models express that upon stretching, the biopolymers of the cytoskeleton – a sub-membranous network of spectrin connections for RBCs – first uncoil, providing relatively little resistance, but when completely stretched out, become practically non-extensible. Between two connected nodes and, the FENE attractive potential reads: (21) where is the maximal distance, and the stretching constant. The force derived from this is: (22) FENE springs exert purely attractive forces. In order to account for the (limited) incompressibility of the spectrin, a simple power law is used (power): (23) The incompressibility coefficient can be derived for the assumption that the total force must vanish for, the equilibrium distance: (24) In the present model, we set, as suggested by [3]. It is convenient to denote the maximal stretch by, the fraction of maximal extension and equilibrium distance. In addition to this purely elastic potential, we also include dissipation as per the Kelvin-Voigt model by adding a parallel dashpot with the damping constant: (25) Here, is the projection of the relative velocity of a pair of connected cortex nodes on their connecting axis. The force is also applied in the direction of the connection. Whereas in-plane stretching and compressive forces can be calculated purely based on the distance between two neighboring cortex nodes, bending forces are calculated for two neighboring triangles. The bending moment between two adjacent triangles is given as (26) Here, is the model parameter determining the bending rigidity, is the instantaneous angle and the spontaneous angle between a pair of triangles with a common edge. A corresponding force is applied to the non-common points of each of the two triangles, with a compensating force applied to the points on the common edge, ensuring that the total force on the cell remains unchanged. This type of bending-stiffness is commonly found in the literature for RBC models, eg. by [4] and [36] - a more general analysis is provided by [37]. Additionally, both a global and local area constraint is used, making sure that both the individual triangle areas and the total area of the red blood cell cannot strongly increase or decrease. As described by [4], this is achieved by a local force with magnitude: (27) in which is the triangle area, the resting triangle area and the local constraint constant. The magnitude of the global force is formulated as: (28) where is the total RBC area, the total resting area and the global constraint constant. For both constants, values were taken from [3]. These forces are applied in the plane of each triangle in the direction from the barycenter of the triangle. Finally, we add a volume constraint since for short timescales, the total cytosol volume of the cell can be considered constant. As for the area, magnitude of the force takes the form (29) with the instantaneous cell volume and the initial cell volume. This force is applied to each node of the cell in its outward direction as found by the Laplace-Beltrami operator, see equation 14. In the low Reynolds number environment in which cells live, motion is dominated by viscous forces [38]. In other words, inertial forces are negligible. For each integration node, Newton' s second law (with explicit Stokes' drag) (30) by leaving out the inertial term, becomes (31) The total force on node is the sum of all the individual forces: Firstly, the forces that are calculated on the triangles are transferred to the nodes – the contact forces only exist for triangles, which are in contact with the substrate. Also, the local and global area constraints for the membrane are added here. Secondly, the cortex connection forces between node and all fixed connections are added, and finally the volume constraint and the gravitational force as well as a random force The total force on node i is the sum of all the individual forces: Firstly, the forces that are calculated on the triangles are transferred to the nodes – the contact forces Fcontact only exist for triangles, which are in contact with the substrate. Also, the local and global area constraints for the membrane are added here. Secondly, the cortex connection forces between node i and all fixed connections k are added, and finally the volume constraint and the gravitational force Fcontact as well as a random force Fcontact for taking into account fluctuations of the membrane can be added. Since those fluctuations do not much influence the spreading dynamics in our simulations, we neglect that term for the results presented. For the right-hand side, we not only discard the term proportional to mass, but we also more explicitly state the components of the constant: starting with the dissipative/friction term generated from the encompassing sphere - substrate friction between two contacting triangles. This coefficient is weighted by the distance of the node from the contact point in that triangle. This ensures symmetry of the friction-matrix (see below) and corresponds to the distribution of the contact force. The component of the substrate friction for a triangle is defined as (compare to e. g. [39]) where is the area of contact in that triangle, is the normalized direction vector between the two encompassing spheres and are, respectively, the normal and tangential friction constants. Secondly, we have the dissipative dashpot of the connections of this node, and lastly we add the drag coefficient for the whole cell in plasma: here, in a first order approximation, we simply divide the formula from Stokes' law by the number of nodes per cell, thereby recapturing the exact result for a spherical cell in Stokes flow. For nodes, whose surrounding triangles are all in contact with the substrate, we define a very high friction constant, effectively fixing those nodes in place. We found that this has no influence on the spreading curves (it can be completely left out), but helps to dampen out small numerical fluctuations in the stiff potential of the contacting plane. This allows us to use larger time steps when solving this equation of motion. Equation 31, which is used in essentially the same form by e. g. [13], [39]–[43], is a first order differential equation, which couples the movements of all particles together. When writing the whole system as (32) it can be shown [13], that the matrix is positive definite, and therefore we are able to solve the system iteratively for the velocities by using the conjugate gradient method. Subsequently, the nodes' movement is integrated by a forward Euler scheme [44]. For a low Reynolds number environment, the amount of kinetic energy (or motion) directly corresponds to the amount of dissipated energy. Equation 32 shows all dissipative terms in the matrix dictating the degree of motion induced by the forces. Identifying all significant dissipative mechanisms is therefore crucial for calculating the dynamics of this system.
How cells spread on a newly encountered surface is an important issue, since it hints at how cells interact physically with the specific material in general. It has been shown before that many cell types have very similar early spreading behavior. This observation has been linked to the mechanical nature of the phenomenon, during which a cell cannot yet react by changing its structure and behavior. Understanding in detail how this passive spreading occurs, and what clues a cell may later respond to is the goal of this work. At the same time, the model we develop here should be very valuable for more complex situations of interacting cells, since it is able to reproduce the purely mechanical response in detail. We find that spreading is limited mainly by energy dissipation upon contact and later dissipation in the cell' s cortex and that no irreversible deformation occurs during the spreading of red blood cells on an adhesive surface.
Analysis of Initial Cell Spreading Using Mechanistic Contact Formulations for a Deformable Cell Model