A viscoelastic–stochastic model of the effects of cytoskeleton remodelling on cell adhesion

Cells can adapt their mechanical properties through cytoskeleton remodelling in response to external stimuli when the cells adhere to the extracellular matrix (ECM). Many studies have investigated the effects of cell and ECM elasticity on cell adhesion. However, experiments determined that cells are viscoelastic and exhibiting stress relaxation, and the mechanism behind the effect of cellular viscoelasticity on the cell adhesion behaviour remains unclear. Therefore, we propose a theoretical model of a cluster of ligand–receptor bonds between two dissimilar viscoelastic media subjected to an applied tensile load. In this model, the distribution of interfacial traction is assumed to follow classical continuum viscoelastic equations, whereas the rupture and rebinding of individual molecular bonds are governed by stochastic equations. On the basis of this model, we determined that viscosity can significantly increase the lifetime, stability and dynamic strength of the adhesion cluster of molecular bonds, because deformation relaxation attributed to the viscoelastic property can increase the rebinding probability of each open bond and reduce the stress concentration in the adhesion area.


Introduction
As a salient characteristic in biological systems, cell adhesion plays a pivotal role in cell differentiation, migration and growth [1][2][3][4][5][6][7]. Moreover, cell adhesion is crucial in the development and maintenance of tissues [8]. Changes in cell adhesion are associated with many diseases [8,9], such as arthritis [10], cancer [11], osteoporosis [12] and atherosclerosis [13]. Generally, most animal cells survive by adhering to the extracellular matrix (ECM) or other cells [14]. Once the cells have detached from the adhesive substrate, apoptosis may be triggered with the dissociation of receptor-ligand bonds, eventually force, implying that whether or not the viscoelasticity affects bond lifetime seems to depend on the bond number [64,65]. However, this conclusion lacks sufficient biophysical supports.
In this paper, we attempt to address the above problem by proposing a coupled viscoelastic-stochastic model of the cell with viscoelastic properties adhering to an elastic substrate under an external applied force, considering that the mechanism of the effect of cellular viscoelasticity on adhesion behaviour remains unclear. On the basis of this model, we studied the coupled effect of cellular deformation, stress/displacement relaxation, stochastically reforming and breaking of receptor-ligand bonds on adhesion strength and stability.

Model
To illustrate the viscoelastic-stochastic adhesion model in this study, we consider the plane strain problem of a linear elastic half space with Young's modulus E s and a linear viscoelastic Kelvin-Voigt half space with Young's modulus E c and viscosity coefficient η c linked by a cluster of ligand-receptor bonds under a remotely applied force F, as shown in figure 1. In addition, the x-axis is placed along the interface between the two half spaces, and the receptor-ligand bonds are uniformly distributed along the x-axis, forming an adhesion patch of size 2a with bond spacing b along both the x-axis and the outof-plane directions, which corresponds to a number density of bonds ρ LR = 1/b 2 . Consider the closed receptor-ligand bond as a linear spring with stiffness k LR , rest length l b and reacting radius l bind around the binding site.

Stress relaxation-dependent interfacial traction
To describe the stress distribution within the adhesion area, we consider the number density of closed bonds at position x, and time t as ρ LR ρ(x, t), where ρ(x, t) is the normalized bond density. Thus, the interfacial stress can be given by where ξ (x, t) is the elastic extension of the bond at position x and time t. If we define the relative normal surface deformation as the sum of surface displacements of two half spaces at the closed-bond position x, w(x, t) = w c (x, t) + w s (x, t), as shown in figure 2, we should have w(x, t) + ξ (x, t) + l b + l bind = W(t), (2.2) where subscripts 'c' and 's' denote cell and substrate, respectively, and W is the interfacial separation between the opposing half spaces at infinity. If the bond at position x and time t is open, then as shown in figure 2, ϕ(x, t) = W(t) − w(x, t) becomes the surface separation at this location. According to theory of continuum mechanics, the relative normal surface deformation can be related to the interfacial stress as (see appendix A for details) is the Laplace transform of the normalized combined creep compliance,J(t) = (E s + E c )J(t). According to Qian et al. [45], α in equation (2.5) can be identified as a dynamic stress concentration index to control the stress-relaxation-dependent distribution of interfacial stress. As the integraldifferential equation ( Figure 2. Schematic of the geometrical relationships at position of (a) closed bond and (b) breaking bond. W(t) is the surface separation at infinity.
theory, thus, similar to the adhesion system of elastic bodies, for the adhesion of viscoelastic medium, a transition from uniform distribution to singular distribution also occurs on the interfacial traction in the adhesion domain when the factor α changes from zero to infinity. For example, when α → 0, from equation (2.4) we can derive w(x, t) ≈ f (t), which means that w(x, t) is a constant at any time t, not a function of spatial coordinates. If we further assume that the bonds are uniformly distributed, then the interfacial traction also becomes constant, corresponding to the state of equally shared loading. However, when α → ∞, solution of equation (2.4) implies γ (x, t) ∝ 1/ 1 − x 2 /a 2 , which is singular at two edges of the adhesion patch at any time t.
On the other hand, considering limiting cases in terms of viscosity, α in equation (2.5) can be rewritten as α 0 = aρ LR k LR [3/(2E s ) + 1/E c ], for η c → 0, and α ∞ = 3aρ LR k LR /(2E s ), for η c → ∞, which are linearly proportional to the adhesion size, bond stiffness and density, and inversely proportional to the combined elastic modulus of the cell and substrate. The existence of α ∞ < α 0 indicates that increased viscosity can reduce the stress concentration in the adhesion domain. If we consider the ratio α 0 /α ∞ = 1 + 2E s /(3E c ), we can find that viscosity can reduce stress concentration index E s /E c times. Even if Young's modulus of the substrate becomes similar to the cell, there is still α 0 /α ∞ ∼ 1.67. This fact indicates that both the cellular elasticity and viscosity are crucial features used by cells to mediate their biological behaviour.

Stochastic dynamics of single receptor-ligand bonds
A closed single receptor-ligand bond shows a finite lifetime only, which inevitably undergoes a transition from the initial closed state to an open state as a result of thermally activated bond dissociation, even in the absence of an external force. Particularly, the dissociation rate of a closed bond located at x n is assumed to increase exponentially with the force F n acting on the bond as [33] where k 0 is the spontaneous dissociation rate in the absence of force and F b is a force scale typically in the pN range; 1/k 0 is suggested to be in the range from a fraction of a second to 100s for receptor-ligand bonds in cell adhesion [66,67]. Additionally, the applied force on the closed bond at position x n can be obtained by F n = k LR ξ (x n ,t).
For free receptor and ligand molecules, bond reform would occur at binding site if the receptor comes sufficiently close to its complementary ligand to react. The bond association rate is assumed to decrease with the separation ϕ between two half spaces as [45,68,69] where k 0 on is the spontaneous reaction rate between the receptor and ligand, k B is the Boltzmann constant, T is the absolute temperature, l bind is the reacting radius around the binding site and Z is the partition function for a receptor confined in a truncated harmonic potential as shown in appendix B.
It can be seen from equations (2.6) and (2.7) that the bond reaction rates are governed by the forces acting on the closed bond and the surface separation for the open bond, respectively. These force and surface separation can be determined through theory of continuum mechanics as shown in appendix C.

Stochastic-viscoelastic coupling through Monte Carlo simulation
By coupling the deformation of the cell and substrate and stochastic dynamics of molecular bonds, we investigate the spatial and temporal evolutions of receptor-ligand bonds using the 'first reaction method' derived from the Gillespie algorithm [70,71]. The key idea of such a method is to sample random trajectories of the system. Specifically, for the adhesion system as shown in figure 1 with total N bonds, the detailed simulation process is described as follows: (1) At time step 0, i.e.t = 0 wheret = k 0 t is the normalized time, we set all N bonds closed. For a given force F, unknown acting forces F m at each binding site x m , m = 1, 2, . . . ,N, can be obtained based on the force balance relation and the associated force-deformation expression for viscoelastic media, as derived in appendix C. (2) Calculate the reaction times at individual bond locations x m by [70,71]

Results and discussions
In this study, bond spacing b is 32 nm [45], stiffness of each receptor-ligand bond k LR is 0.25 pN nm −1 [68], rest length of each receptor-ligand bond is 25 nm [56,72], force scale F b is 4 pN [38,39], and a 0 = 5 nm denotes the radius of individual bonds [73]. The spontaneous dissociation and association rates are taken    as k 0 = 5 s −1 [43] and k 0 on = 500 s −1 [66], respectively. The reaction radius of the receptors and ligands is set to 1 nm [45]. Young's modulus of the substrate and cell is 150 kPa and 10 kPa, respectively.
The relevant parameters adopted in this study are summarized in table 1. For the former, in the light of the spatial and temporal averages, figure 4b presents the steady-state distribution of the interfacial stress for the adhesion of viscoelastic bodies with different viscosities through a cluster of molecular bonds and under a remotely applied load. Figure 4b indicates that larger viscosity corresponds to a more uniformly distributed interfacial stress. From the view of individual bonds, large loading force would result in unstable binding between receptor and ligand molecules based on Bell's theory [33]. As shown in figure 4b, the discrete interfacial force in adhesion edges is larger than that in adhesion centre. Therefore, there is greater probability that the bonds break at edges. Once the bonds at edges break, the crack would propagate from edge to centre. On the contrary, if the force has a uniform distribution, all bonds would have stable state for adhesion.   For the latter, we note that there exist characteristic time scales for the behaviours of viscoelastic relaxation, bond dissociation and association, which can be represented by T 1 = η c /E c , 1/k 0 and 1/k 0 on , respectively. Competitions between these time scales determine how the cellular viscoelastic property may affect the adhesion behaviour of individual ligand-receptor bonds. For a closed bond, if the characteristic bond dissociation time 1/k 0 is longer than that of the viscoelastic relaxation time T 1 , the deformation relaxation has no influence on bond breaking. For an open bond, when the characteristic bond association time 1/k 0 on is longer than the characteristic relaxation time T 1 , the bond association is also independent of the cell viscosity, but when 1/k 0 on < T 1 , stress relaxation may postpone the recovery of interfacial deformation from small surface separation, which makes the reforming event easier to occur. In typical biological system, Young's modulus of cells, E c , is of the order of 10 kPa or less [74], the cell viscosity ranges from hundreds to thousands of Pa s [75], the spontaneous bond dissociation rate is in the range of 5 × 10 −6 -9 s −1 [76,77] and the spontaneous bond association rate is in the range of 1-10 3 s −1 [76]. Therefore, in most cases there are always 1/k 0 > T 1 and 1/k 0 on < T 1 , which implies that bond breaking is usually independent of cell viscosity, but the bond rebinding is enhanced by the cell viscosity.

Influence of cellular viscosity on the lifetime of receptor-ligand bond cluster
In order to provide an intuitive impression of how the viscosity may increase the bond rebinding possibility, we consider an illustrative example that there is a single bond which links the cell and substrate, as shown in figure 5. If the cell has high viscosity, the recovery procedure for the cell deformation must be slower than that of the cell deformation with low viscosity after bond breaking, making a small interfacial separation between the cell and substrate. This observation indicates that the next reaction of bond rebinding at the reaction site for a cell with high viscosity occurs easier than that for the cell with low viscosity. In an extreme case, if the cell is purely elastic, deformation of the cell can be immediately recovered as long as the bond breaks. Hence, for the elastic system, large interfacial separation associates with low probability of bond rebinding.
From a biophysical point of view, cells may take advantage of their own mechanical properties for proper functioning.    Figure 6 shows that lifetime strongly depends on cell viscosity with small applied forces although this dependence becomes weak at large applied forces. This phenomenon occurs because if the applied force is large, the interfacial separation also increases. Thus, the rebinding events will rarely occur. Gupta [65] considered the adhesion of a polymorphonuclear leucocyte cell. He used a spherical membrane containing a Newtonian fluid to model the cell, and the cell was adhered to the substrate by a single molecular bond. Based on this model, Gupta [65] concluded that cellular viscoelasticity does not affect the average bond lifetime for constant-force loading, because the bond force becomes equal to the applied constant force after a transient period. This conclusion is unfortunately biased. As shown in figure 6, the rupture events of the bonds dominate and the rebinding events become rare only if the loading force is very large, and then the bond lifetime is independent of cellular viscosity. When the loading force is relatively small, the cellular viscosity will significantly influence the bond lifetime. In the work by Gupta [65], the loading force is of the order of hundreds of pNs, which can be more than 30 times larger than the typical force scale F b , leading to a very large rupture rate according to equation (2.6). By contrast, we chose the loading force in tens of PNs, which can lead to a more reasonable rupture rate. Under the present force range, the bond lifetime strongly depends on cellular viscosity, as indicated in figure 6.

equation (2.3) can be reduced to
The above equation implies that in the case of η c → ∞, only the substrate would deform in response to the applied force acting on the interface. For η c → 0, both the substrate and cell would deform. Previous experimental and theoretical studies [23,24,45] demonstrated that substrate elasticity can regulate cell adhesion. In this study, we broaden this conclusion by determining that the cellular viscosity can also significantly affect cell adhesion.

Effects of cellular viscosity on the dynamic strength of bond cluster
Our model can also be expanded to predict the dynamic adhesion between the cell and substrate via receptor-ligand bonds. We consider the displacement controlled loading by setting interfacial separation h as h = λt, where λ is the loading rate. In this case, consultant interfacial force F is no longer constant. The mean rupture force of a bond cluster with N = 50 is numerically determined and plotted in figure 7 as a function of the loading rate under different cell viscosities. The figure indicates that, if the loading rate is low, the mean rupture force exhibits an asymptotical strength limit and the cellular viscosity only slightly affects the bond strength. This phenomenon was also observed by Li & Ji [78,79] in studying the stretch of a single molecular bond. The asymptotical bond strength results from the reaction equilibrium between the receptor and ligand, as most recently revealed by Li et al. [80]. Furthermore, if the loading rate is large, the dynamic strength of the bond cluster increases with loading rate and the cellular viscosity slightly affects adhesion strength.
If the loading rate is intermediate large, as shown in figure 7, cellular viscosity can significantly enhance the dynamic strength of the bond cluster. This loading rate range that leads to viscositydependent bond strength is a result of the competition between time scales of bond formation, loading rate and creep. When the cell separates from the substrate under a low loading rate, interfacial deformation has sufficient time to recover so that the cellular viscosity does not affect mean rupture force. In the case of large loading rate, the dynamic strength is dominated by the fast increase of interfacial separation between the cell and substrate, implying that the rebinding events rarely occur after bondbreaking. Thus, the rupture force of the bond cluster depends on the cellular viscosity only in the intermediate range of loading rate.
Moreover, figure 7 shows the upper (η c → ∞) and lower (η c → 0) bounds on the curves of effective mean rupture force versus loading rate for different viscosities. The figure indicates that the rupture force of the bond cluster does not change with viscosity when the viscosity exceeds 50 kPa s. Interestingly, in nature, the value of viscosity for most cells is in the range of tens to thousands of Pa s [81][82][83], implying that cells may adapt their own viscosity by cytoskeleton remodelling to control adhesion.

Effects of cellular viscosity on the window of stable adhesion
Normalized lifetime of the adhesion cluster is shown in figure 8 as a function of cluster size 2a/b for constant stress F/2ab = 0.53 kPa and different cell viscosities. Evidently, a size-window exists for the stable adhesion of the molecular bond cluster. These size effects on cell adhesion were demonstrated based on the model that disregards the effects of viscosity [40,45,50]. By contrast, in this study, we found that cell viscosity can increase the lifetime and also broaden the size of the stable window, indicating that the cell may use its own viscous property other than or together with elasticity property to maintain the stable adhesion in a large range of cluster sizes. Why there exists such a size-window: if bond cluster size is small, since few bonds are involved in maintaining stable adhesion, the stochastic breaking and reforming of molecular bonds would lead to a short lifetime. For very large bond cluster size, the stress concentration results in the cluster breaking from edge to centre and also eventually makes a short lifetime. Therefore, there will be an optimal cluster size that maximizes the adhesion strength.

Conclusion
In this research, we developed an idealized viscoelastic-stochastic model of two elastic and viscoelastic bodies connected by a receptor-ligand bond cluster. Based on this model, we have performed Monte Carlo simulations to determine the role of cellular viscosity in biological behaviour of cell adhesion, such as lifetime, dynamic strength and stable window, by coupling the continuum deformation of the viscoelastic cell and elastic substrate and the stochastic behaviour of molecular bonds. The main conclusions are summarized as follows. (1) We broadened the concept of the stress concentration index given by Qian et al. [45] to include the effect of stress relaxation. Based on such a dynamic stress concentration index, we pointed out that cellular viscosity can reduce stress concentration at the adhesion domain. (2) We sufficiently demonstrated that the time-dependent recovery process of viscoelastic deformation increases the probability of rebinding for free receptors and ligands at binding sites. Thus, cellular viscosity can generally prolong the lifetime of an adhesion cluster. However, once the loading force is large, the effect of cellular viscosity on the bond lifetime becomes negligible. That is why Gupta [65] did not observe the effect of cellular viscosity on bond lifetime under constant-force loading.
(3) Compared with the results in the previous study [64], our analysis clearly shows that cellular viscosity enhances rather than reduces the adhesion strength of a cluster of molecular bonds under a displacement controlled dynamic loading. (4) Cellular viscosity also plays an important role in controlling the sizewindow for stable adhesion of a molecular bond cluster. Cellular viscosity can effectively broaden the stable size-window and prolong cluster lifetime but only for clusters within this size-window. This finding indicates that cellular viscosity can exert the same effect on the adhesion stabilization as substrate stiffness [45].
In vitro experiments on cells can be designed to verify the proposed theory, for example, some nanoparticles and drugs can be used to change the viscosity of the cells, and such a viscosity change can be reflected through the change in their adhesion states as revealed by the present theory.
Overall, the results based on this model, although overly simplified in this study, seem to support a novel view on understanding the role of cell viscosity in cell adhesion. For such linear problem, in Laplace-transform space, the surface displacement gradient ∂ŵ c (x, s)/∂x of the viscoelastic body can be related to the interfacial tractionγ (x, s) through Green's function where a is half size of the adhesion patch,γ (x, s) andŵ(x, s) correspond to the Laplace transform of γ (x,t) and w(x, t) with respect to time t. Applying inverse Laplace transform to equation (A 13) yields Meanwhile, the deformation of the elastic substrate would be induced by the interfacial traction γ (x,t). The surface displacement gradient ∂w s (x, t)/∂x of the elastic body has a relation with the interfacial traction γ (x, t) through Green's function as [86] where the relative normal surface deformation w(x, t) is denoted as the sum of surface deformations of the elastic and viscoelastic half spaces at position x and time t, i.e. w(x, t) = w c (x, t) + w s (x, t), as shown in figure 2. If the cell is considered as the Kelvin-Voigt solid with Young's modulus E c and viscosity coefficient η c , then the combined creep compliance J(t) can be expressed as where J s (t) = 3/2E s , J c (t) = (1 − exp(−t/T 1 ))/E c and T 1 = η c /E c represent the characteristic time scale of relaxation.

Appendix B. Association rate of receptor-ligand bonds
The reaction between receptor and ligand molecules is mediated by the surface separation ϕ as shown in figure 2. At a given separation ϕ, the truncated harmonic potential for a free receptor molecules has the form where y is the possible displacement of the receptor molecule as displayed in figure 2. Consider that the receptor displacement follows the Boltzmann distribution, the probability density function for the receptor with displacement y is where Z is the partition function, which can be given by solving the normalization equation where erf(.) is the error function. Hence, the probability that the receptor falls into a reacting radius l bind around the binding site is where k 0 on is the spontaneous reaction rate between the receptor and ligand.

Appendix C. Forces and displacements for discrete bonds
For the cell-substrate adhesion system as shown in figure 1, we consider the stochastic process of the cluster of discretely distributed receptor-ligand bonds. We use theory of continuum mechanics to determine the forces and displacements at bond sites. Specifically, for a given bond location x m within the adhesion patch, we will determine relative interfacial displacement w mn induced by a force F n acting on another bond location, x n (m = n).
In appendix A, we have derived Integrating both sides of equation (C 1) yields where we have assumed w(x ∞ ,t) = 0 when x ∞ is sufficiently large. Assuming that only one bond force, F n (t), is applied at bond location x n , so that On the other hand, in order to derive the displacement at position x n to be induced by the bond force F n at the same position, we replace the concentrate force F n by an equivalent uniform pressure with half-width a 0 . In this case, equation (C 2) gives w(x n , t) = 1 2a 0 bπ and x ∞ x n , x ∞ a 0 have been considered. Therefore, the self-displacement at x n induced by the force F n can be given by Then, for a given time t, the total relative interfacial displacement at location x m becomes w(x m ) = w mm + N n=1,n =m w mn where N is the total number of receptor-ligand bonds. For open-bond location x n , as shown in figure 2 the separation between two half spaces, ϕ, at position, x n , presents the form ϕ(x n , t) = W(t) − w(x n , t), (C 11) and the exerted force is zero, F n = 0. For closed-bond location x m , the geometrical relationship is w(x m , t) + ξ (x m , t) + l b + l bind = W(t).
(C 12) Combining equations (C 10)-(C 12) gives the exerting force on each bond and the separations for open bonds.