{'fraction_non_alphanumeric': 0.025983980377879124, 'fraction_numerical': 0.013375234737285862, 'mean_word_length': 4.967638650657519, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling and simulation of chemomechanics at the cell-matrix interface

Chemomechanical characteristics of the extracellular materials with which cells interact can have a profound impact on cell adhesion and migration. To understand and modulate such complex multiscale processes, a detailed understanding of the feedback between a cell and the adjacent microenvironment is crucial. Here, we use computational modeling and simulation to examine the cell-matrix interaction at both the molecular and continuum lengthscales. Using steered molecular dynamics, we consider how extracellular matrix (ECM) stiffness and extracellular pH influence the interaction between cell surface adhesion receptors and extracellular matrix ligands, and we predict potential consequences for focal adhesion formation and dissolution. Using continuum-level finite element simulations and analytical methods to model cell-induced ECM deformation as a function of ECM stiffness and thickness, we consider the implications toward design of synthetic substrata for cell biology experiments that intend to decouple chemical and mechanical cues.


Introduction and Background
Chemomechanics of the extracellular matrix. Critical cell processes such as matrix adhesion, proliferation, migration, differentiation and intracellular signaling are strongly influenced by chemomechanical properties of the surrounding cellular microenvironment. For example, cell functions can be altered by chemical factors such as local pH changes or the addition of soluble growth factors and synthetic drugs. [1][2][3][4] Cells also respond to mechanical stimuli defined by or transmitted through the extracellular matrix (ECM) and convert this input into chemical signals. The detailed mechanisms of this chemomechanical feedback are unclear and elucidation of the key structures and cue-signal-response networks is currently an active area of research.
Mechanical forces at the cell-matrix interface can be transmitted through protein complexes called focal complexes (FCs) or focal adhesions (FAs) (Fig. 1A). Focal complexes are considered small (<1 μm 2 ) dot-like structures that are commonly involved in cell migration. 5,6 Focal adhesions are more spatiotemporally stable and include more protein components, and are thus larger in size (≥1 μm 2 ). 5 These macromolecular complexes form the physical linkage between the cell and the ECM through individual interactions between an integrin receptor and an ECM ligand 6 (Fig. 1B). Integrins are heterodimeric transmembrane receptors that can bind to several distinct protein ligands comprising the ECM, such as fibronectin, collagen and laminin. 7 On the cytoplasmic side of the cell surface, integrins are linked to the actin cytoskeleton through a multi-protein assembly. 7 As a result of this physical engagement of the cytoskeleton, cells can transmit actomyosin-generated forces against the ECM.
It is known that ECM stiffness spans a considerable range in healthy [8][9][10] and diseased 11 tissues. This is often quantified by the Young's modulus of the ECM material E, the ratio of macroscopic tensile stress to strain in the region of linear elastic deformation. For example, brain tissue is reported to exhibit E of ~2.5 kPa, 10 whereas E of muscle tissue is ~12 kPa 8 and E of trabecular bone tissue is ~18 GPa. 12 Normal breast tissue exhibits E of ~0.15 kPa, but breast tumors can exhibit E as large as ~4 kPa. 11 Although E is a gross approximation of the mechanical response of (bio)polymeric ECM and tissues (which are viscoelastic and, often, highly nonlinear elastic materials), 13 it is increasingly apparent that changes in cell morphology and certain functions correlate with changes in matrix stiffness. 9,11,[14][15][16] Experiments that probe cell response to matrix mechanics often utilize synthetic substrata of controlled mechanical stiffness (or, inversely, compliance). Such materials include polyelectrolyte multilayers (PEMs), 14,17 polyacrylamide hydrogels 18 and protein hydrogels, 19,20 including complex commercial compositions such as Matrigel ® .
Several recent studies have demonstrated that matrix stiffness can affect a wide variety of cell processes. Using polyelectrolyte multi-layers, Thompson et al., showed that human dermal microvascular endothelial cells displayed increased frequency of adhesion on stiffer PEM substrata, over the range of 0.5 MPa < E < 100 MPa. 14 Paszek et al., demonstrated that substrata stiffness can also influence cell signaling. The authors used fibronectin-modified polyacrylamide gels to show that cells which were adhered to stiffer gels exhibited increased activity of the protein kinase ERK and the GTPase Rho, which promotes FA formation. 11 Engler et al., demonstrated that collagen-overcoated polyacrylamide gels with stiffnesses comparable to certain tissues can direct adult mesenchymal stem cells to display some markers for those tissue cell types. 9 In these three examples, cells were cultured in two-dimensional (2D) constructs, i.e., on top of the substrata. This is considered an appropriate in vitro tissue model for epithelial-type cells, which exist in an approximately 2D environment in vivo.
The quantification of stiffness becomes more complicated for a fully 3D matrix with large pores or mesh size, in which a FC/FA might only span one individual component strut of the matrix (Fig.  2). Some materials designed for 3D cell culture fall into this category, such as collagen-GAG matrices used for studies of fibroblast-mediated ECM contraction. 21 With such materials, one could consider the stiffness of an individual component strut or the composite stiffness of the matrix (poroelasticity), which incorporates the effects of the matrix architecture such as %-porosity. It is unclear as yet which measure is more relevant for specific cell responses. Moreover, the effective structural stiffness of an individual strut or pore wall depends on the length and the network continuity of the strut. In one example using 3D cell culture, Zaman et al., demonstrated that DU-145 human prostate carcinoma cell migration speed in Matrigel ® is biphasic with respect to composite matrix stiffness, with maximum cell speed at an intermediate stiffness. 16 However, the ligand concentration in the matrix also varied with stiffness, and the authors noted that maximum cell speed depended on the optimum balance between these chemical and mechanical parameters.
The above studies demonstrate a variety of synthetic substrata/ matrices and measured indicators of cell function, but represent only a small sample of the published work on substrata mechanics and its influence on cell behavior. With any study of this nature, it is also important to understand that the interplay between chemical and mechanical feedback is critical in determining cell response in a physiological environment.  Cell adhered to a 2D substratum of pore interstitial spacing smaller than cell width. A focal adhesion or focal complex can span many struts and the composite material stiffness may be an appropriate characterization; (B) Cell embedded in a 3D matrix with pore interstitial spacing greater than cell width. A focal adhesion or focal complex can only contact one strut and proper characterization of stiffness depends on strut material stiffness, length and network connectivity. experiments, it is possible to use rupture force (F R ) values obtained at different loading rates to extract equilibrium kinetic and energetic parameters of unbinding through the application of modified Bell's models. 23,24 However, in order for the simulation to be computationally feasible with current computational resources, SMD loading rates significantly exceed experimental loading rates. In such cases, equilibrium kinetic parameters obtained through application of Bell's model will not quantitatively match experimental parameters. 25,26 However, it is possible to construct an equilibrium free energy profile from nonequilibrium SMD simulations through the application of Jarzynski's equality. 27 Furthermore, SMD is useful for measuring relative values of rupture forces or times, and this technique allows quantitative study of unbinding mechanisms and of the hierarchy of weak and strong interactions. In addition, SMD enables us to study in detail the range of responses that can result from different initial conditions and/or molecular configurations.
SMD has been used to study protein unfolding as well as ligand-receptor unbinding. One of the most notable unfolding studies concerns the muscle protein titin. A SMD study revealed the forcebearing region of the titin structure, 28 and this finding was supported by an experimental study of a titin protein mutant that weakened this interaction. 29 SMD has also been very useful in the study of the extracellular matrix protein fibronectin. Here, simulations have shown that applied force can regulate accessibility of the Arg-Gly-Asp (RGD) binding site 30 and of the synergy site in fibronectin that enhances binding to some integrins. 31 Such simulations have been used to probe the hierarchy of mechanical stability among several fibronectin domains, 32 and predictions have been supported by atomic force microscope (AFM) molecular force spectroscopy (MFS) experiments. 33 Furthermore, SMD has indicated the existence of multiple unfolding pathways for the fibronectin FN-III 10 and FN-III 1 domains, [34][35][36] and AFM MFS results are generally in agreement with these findings. 36,37 In the context of ligand-receptor unbinding, SMD has long been used to simulate the biotin-streptavidin complex, a molecular pair of high and specific binding affinity, in order to calculate rupture forces and examine possible unbinding mechanisms. 22,24 Improvements in computational resources have since allowed for longer and more accurate simulations of this system. 24,38 SMD has also been used to study the unbinding of the hormone retinoic acid from its receptor, and results suggest that multiple un/binding pathways exist. 39 SMD has also elucidated mechanisms for the binding of retinal to bacterio-opsin, which is an important step in the formation of bacteriorhodopsin. 40 In addition, fluorescein-antifluorescein antibody unbinding was characterized with SMD, and it was shown that the unbinding mechanism changed upon certain mutations of the antibody. 41,42 For our studies of ligand-receptor chemomechanics, we first consider the biotin-streptavidin complex to test the effect of initial bound configuration and ECM stiffness on unbinding. Biotinstreptavidin is a ligand-receptor system that has been very well characterized both experimentally and computationally, 22,38,[43][44][45][46][47][48][49][50][51] making this complex an ideal model system for fundamental ligandreceptor studies. Here, we use the biotin-streptavidin results to inform our design of a study for the integrin-RGD complex that is more relevant to cell-substrata adhesion.
In this paper, we focus on the following central questions: How do chemomechanical properties of the extracellular environment influence the cell-matrix interaction; and how does this feedback affect cell adhesion, substratum deformation and migration potential? We concentrate mainly on the effect of substrata stiffness and apply this question at multiple length scales of the cell and of the substratum. At the atomistic/molecular scale, the unit physical linkage between a cell and the ECM is a ligand-receptor interaction, and the energetics and kinetics of this reversible binding likely play an important part in the chemomechanical feedback often termed mechanotransduction. At this scale, we seek to understand how matrix stiffness affects the integrin receptor-ECM ligand interaction, and in particular the binding lifetime of the key integrin-RGD complex. Such mechanistic understanding should elucidate how ECM stiffness affects the ability of FCs or FAs to form, stabilize and grow. At the continuum scale, we seek to understand how a cell transmits forces to and deforms the adjacent substrata after mature FAs have formed. In other words, over what lengthscales does a cell probe and sample its mechanical microenvironment? Effective design of synthetic substrata that accurately mimic the physiological niche requires this continuum level understanding of how cells interact with adjacent materials. This predictive capability should be useful in the consideration of appropriate substrata thickness for cell studies such as those described above. Furthermore, mechanical deformation of the cell and of the substrata/matrix can influence concentration gradients of important chemical species (e.g., growth factors) within the ECM, which could induce localization or sequestration of molecules and thus affect cell signaling and downstream behaviors.
Computational models. Computational methods provide us with an efficient (albeit simplified) means to explore these questions at multiple lengthscales. Experimentally, it can be challenging to manipulate mechanical properties such as substrata stiffness without concurrently changing other important system characteristics, but computational models can decouple such variables systematically. In addition, the use of simplified model systems can reduce computational expense but still allow characterization of general principles that can be applied to more specific and complex systems.
In order to study the effects of ECM stiffness on binding properties of ECM ligands and cell surface receptors, we use steered molecular dynamics (SMD) to model unbinding of individual ligand-receptor pairs at the atomistic level. It is not currently computationally feasible to simulate an entire focal adhesion complex with this level of detail, but we can gain important knowledge and form reasonable coarse-graining approaches by studying individual ligandreceptor pairs. SMD enables computational modeling of conditions comparable to single molecule force spectroscopy experiments performed with tools such as the atomic force microscope or optical tweezers. With SMD, an external force is applied in order to displace the ligand from the receptor cleft or binding pocket. 22 Thus, it becomes possible to observe unbinding or unfolding events over the restricted timescales currently available to atomistic simulation techniques such as molecular dynamics (<1 μs). Typically, external force is applied with a force transducer of known stiffness k c , a computational spring that is moved at constant velocity. The force applied to the complex F can be calculated from the spring extension x, e.g., F = k c x for a Hookean spring. With molecular force spectroscopy To study cell-induced deformation of the ECM at the continuum level, we use finite element modeling (FEM) to represent a stable FA-ECM complex. Such computational modeling allows us to easily vary ECM mechanical properties, as well as substrata geometric properties such as thickness, without perturbing the chemistry of the system. FEM is very useful for understanding material deformation in response to force, and has been used to simulate many different biological systems. FEM often employs simplified representation of cell shape and generally represents biological materials as homogeneous continua defined by constitutive laws of linear (visco)elastoplastic deformation. However, this modeling method also allows for more complex representations when appropriate. 52 Such models can be used to gain understanding of cell-level biomechanical experiments, just as SMD lends insight into single-molecule experiments. For example, FEM has been used to understand cell deformation during magnetic bead twisting experiments, 53 and to simulate deformation of a red blood cell with optical tweezers. The latter model allowed the authors to vary factors such as cell size, membrane viscosity and applied force to observe the resulting effects on macroscopic cell deformation. 54 In another study, atomic force microscopy was used to measure mechanical properties and topography of cell components. This information was then incorporated into a finite element model that includes a complex cell shape and heterogeneity of properties. This model was used to determine cellular strains in response to several whole-cell and micromanipulation experimental techniques. 55 In addition to simulating and defining cell mechanical experiments on cells, FEM has been used to represent in vivo mechanical environments. Such models have been used to study the influence of an embedded chondrocyte on stress-strain fields imposed on the surrounding ECM, 56 and to model cell deformation during leukocyte rolling. In the latter study, different cell membrane stiffnesses were considered in order to compare the resulting deformation and rolling behavior. 57 FEM can also be used to study interactions between two cells or between two masses of cells. 58 Finally, although FEM is commonly used to study macroscopic structures such as cells, they can be used to represent coarse-grained proteins as well. In an interesting multiscale study, FEM was used to simulate a mechanosensitive channel in a lipid membrane. Molecular statics calculations were used to determine parameters for interactions between the finite element model components. With this model, the authors were able to study how different types of membrane deformation influence channel opening. 59 In our studies, we use a finite element model of focal adhesions adhered to an ECM in order to study ECM deformation due to FA traction, as a function of substrata stiffness and thickness. Our objective is to understand the lengthscales over which cell-generated traction deforms the underlying substrata, and to observe how deformation is influenced by geometric and mechanical properties of the ECM. We also compare the FEM results to an analytical model drawn from continuum elasticity theory.

Clustering analysis is an effective tool for selection of input configurations for SMD simulations.
It is an accepted practice for SMD studies to use a single initial bound configuration for simulation and to draw general conclusions from the results. 39,40,66 However, as computational power increases, it is becoming apparent that it is desirable to simulate unbinding of multiple initial configurations in order to observe variations in output. 24,41,42 We examined the distribution in rupture force F R of the biotin-streptavidin system by using a large (N > 50) set of equilibrated bound configurations as input for SMD simulations of this biomolecular complex. All of these configurations come from a single equilibration trajectory, and thus we expect that differences in atomic positions among these configurations are subtle. In these simulations, the force transducer had a spring constant of 1686 kJ mol -1 nm -2 (2.8 N/m) and was displaced with a velocity of 0.8 m/s. Simulation results show that the measured rupture force can span a range of ~600 pN, with the average <F R > = 750 ± 139 pN (Fig. 3). A single simulation would give a single value for F R , which is not an adequate characterization of this rupture force distribution.
While simulation of 50+ initial configurations allowed us to examine the rupture force distribution in detail, this number of simulations is very computationally and time intensive. For the purposes of simulating different systems, it is more desirable to characterize unbinding behavior with a minimal number of independent simulations. Clustering analysis is a useful technique for finding conformational subsets and categorizing configurations within a molecular dynamics trajectory, [70][71][72][73][74][75][76][77][78] and we applied this approach in order to identify a small set of input configurations for SMD simulations. We clustered with respect to the RMSD of three separate groups, termed "clustering groups" below: (1) the full ligand-receptor system, including the streptavidin tetramer and all four bound biotin molecules; (2) the residues in one streptavidin binding pocket with its bound biotin; (3) the residues in the streptavidin binding pocket Figure 3. A total of 53 configurations were used as input for SMD simulations (filled circles). The rupture force distribution was compared to results for cluster centroids from three different clustering groups. Group 1: full streptavidin tetramer and four bound biotins (squares); Group 2: streptavidin binding pocket residues and one bound biotin (triangles); Group 3: streptavidin binding pocket residues only (open circles). The cluster centroids from all three clustering groups produce an average rupture force that is not statistically different than that of the full set of initial configurations.
alone. We selected the cluster centroids as representative configurations for input into SMD simulations, and tested whether the resulting rupture forces adequately reproduced the rupture force distribution that we aimed to capture. As shown in Figure 3, the cluster centroids from all three clustering groups result in average <F R > values that are not statistically significantly different from the average <F R > of the full set of 50+ configurations. Therefore, using any of these clustering methods, one can characterize rupture force distribution with a tractable set (N = 10-11) of initial configurations. A rigorous analysis of the properties of the rupture force distribution is ongoing, in order to determine the minimum number of simulations necessary for adequate characterization of ligand-receptor unbinding simulated via SMD (in preparation).
In summary, clustering analysis is an effective method for objectively selecting input configurations for SMD simulations. This approach provides a finite and practical number of configurations to sample, and these configurations adequately characterize the unbinding behavior of the system. As this general method can be applied to any ligand-receptor complex, we use such cluster sampling in our studies of the integrin α v β b -RGD system described below. However, we first examine the effect of force transducer stiffness on biotin-streptavidin binding lifetime.
Biotin-streptavidin binding lifetime varies inversely with tether stiffness. Using Bell's model of ligand-receptor dissociation rates under external mechanical force, one can use F R values obtained at different loading rates F ' to extract the equilibrium binding lifetime. 23 Bell's model implicitly states that at a fixed F ', the lifetime of a ligand-receptor complex should not depend on the force transducer or tether stiffness k c (beyond the role that k c plays in defining F ' = k c v ⇀ . 23 ). In our SMD simulations, k c is the stiffness of a virtual Hookean spring that is used to pull the biotin ligand out of the streptavidin binding pocket. In a biological context, k c could represent the stiffness of the ECM to which a ligand is attached. Counter to this assumption in Bell's model, when we tested k c of 0.83 N/m, 4.15 N/m and 8.30 N/m at a fixed loading rate, our simulation results show that the equilibrium binding lifetime (see Methods and Fig. 4A and B) of the complex τ is strongly dependent on k c and varies inversely with the stiffness of the tether (Fig. 4C). In our simulations, τ varies from 7.60 x 10 6 s for k c = 8.30 N/m to 6.90 x 10 8 s for k c = 0.83 N/m. This result is confirmed by AFM MFS measurements on the biotin-streptavidin complex acquired by our group, 24 and others have recently noted similar, strong effects of k c on F R and inferred energetics for molecular complexes 79 and for self-assembling monolayers. 26 Because the SMD simulations use loading rates in gross excess of the experimentally accessible range, values for binding lifetime do not agree quantitatively with experimental values. However, the trends in behavior match qualitatively (experiments: 5.0 x 10 5 s for k c = 0.058 N/m to 1.20 x 10 6 s for k c = 0.035 N/m 24 ), indicating that the predicted dependence of τ on k c is valid. This dependence is due to the fact that the biotin-streptavidin energy landscape and barrier heights depend on k c . 24 Moreover, simply attaching a harmonic spring to the complex may introduce an additional energy barrier, even with zero applied force. 26 These results prompted ongoing studies in our lab, extending this work to the integrin-RGD system. If the dependence on k c holds true for integrin-RGD, this would suggest that a ligand-receptor complex linked to a stiffer ECM tether would have a shorter lifetime than a complex linked to a more compliant ECM. A change in integrin-RGD binding lifetime could affect the dynamics of FC/FA formation, as it takes a significant amount of time (~60 minutes 5 ) for all of the component proteins to cluster and assemble. This could be important in a pathological context such as a tumor microenvironment, which typically exhibits greater stiffness than normal tissue. 11 Acidic extracellular pH does not significantly alter integrin α v β 3 -RGD rupture force. Our SMD results with the biotin-streptavidin model system have demonstrated that factors such as tether stiffness can have a strong effect on unbinding behavior. With this knowledge, we can apply this methodology to a ligand-receptor system that is essential for cell adhesion and migration, such as integrin α v β 3 -RGD. Here, we are also interested in the chemical environment of the complex, in addition to the mechanical environment. In the environment of a tumor or wound bed, it is well established that the extracellular pH can become acidic. [80][81][82] In the tumor environment, the measured pH values range from 5.85-7.68. 81 Early stage wounds have a measured average pH of 5.7-6.1. 80 This acidic pH could affect the integrin-ligand interaction and influence binding/unbinding energetics and kinetics. This could then alter FC/FA dynamics and cell behaviors that depend on FC/ FA formation/dissolution, such as cell migration.
To simulate an acidic pH with SMD, we used Multi-Conformation Continuum Electrostatics 64,65 to predict accurate pKa values for all the titratable amino acid residues in the system. We protonated all residues with pKa > 6.0 in order to simulate an effective extracellular pH of ~5.5. From the biotin-streptavidin results, we have learned that it is essential to sample many input configurations and observe the range in behavior. Therefore, we chose to generate multiple equilibrated trajectories for each pH, and to use clustering analysis based on RMSD of binding cleft residues to select SMD input configurations from each equilibrated trajectory. In this way, we can sample a wide range of equilibrated configurations. The biotin-streptavidin results also showed that tether stiffness k c can have a strong effect on the measured rupture force. Therefore, we use only a single loading rate F ' (8.28 N/s) and k c (4.14 N/m) value for these initial integrin-RGD simulations. Systematic consideration of independent variations in both F ' and k c are ongoing.
For these fixed F ' and k c values , our SMD results show that the complex at pH 5.5 displays a larger variation in rupture force than at pH 7.4, but that there is no statistically significant difference in average rupture force as a function of pH over this range (Fig. 5). However, there may still be important differences in unbinding mechanisms or among ostensibly equilibrated configurations that are not evident from these rupture force data alone. This is one focus of our current study of this complex.
Multiple interpretations of substrata critical thickness exist. Thus far, our results have provided information about the factors that are involved in ligand-receptor dynamics required of FA formation. Now we shift our focus to understanding cell behavior after mature FAs have formed. More specifically, we would like to understand the length scales over which a cell incurs displacement of the underlying ECM as a function of ECM properties. This has particular relevance to the design of engineered substrata for biological experiments. In these types of experiments, synthetic substrata are designed to represent the relatively compliant mechanical environment of an ECM, but they are often placed on top of comparably rigid materials such as glass or polystyrene. Therefore, if the synthetic substrata are to be an accurate mechanical analogue of the ECM, it is required that the mechanical properties that the cell "feels" via traction-induced strain fields within the ECM are not skewed by the proximity of the underlying rigid material.
Using a finite element model of cell FAs adhered to an ECM, we first examined the von Mises stress distribution for a FA complex adhered to an ECM of thickness 6 μm and E = 1800 kPa or 18 kPa (Fig. 6A and B). Our results show that the stress distribution is clearly a function of ECM stiffness. Stress distribution will affect ECM deformation, which in turn could influence concentration gradients and localization of important chemical factors. This is one explanation for how a change in ECM stiffness could be both correlative and causal in modified cell signaling and resulting behaviors.
We also examined the principal strain field through the depth of the ECM, in response to traction on substrata of E = 1800 kPa and of thickness t = 6 μm, 3 μm or 0.5 μm. Our results show that the principal strain fields decayed to <0.1% through the thickness of this stiff ECM for t = 3 μm or 6 μm. However, for ECM thickness t = 0.5 μm, the strain field did not completely decay within the ECM (Fig.  6C). This indicates that the traction is not completely dissipated by the ECM material, and that the effective mechanical stiffness of the cell substrata includes a composite response with contribution from the underlying rigid support.
For the purposes of engineering substrata, it would be useful to have a measure such as a critical ECM thickness, above which the cell would no longer be influenced by any underlying material. Based on these FEM simulations, one definition of critical thickness is the depth through which the principal strain is ≥0.1%. This value was chosen to represent the point at which the strain has decayed effectively to zero. Using this definition, we examined critical thickness as a function of total thickness and stiffness of the ECM. As shown in Figure 6D, the critical thickness naturally increases as the ECM becomes more compliant. Our results show that a substratum of E = 18 kPa must be at least 3 μm thick for the strain field to decrease to a negligible value near the underlying rigid support. In contrast, a substratum of E = 1800 kPa would only need to be 0.375 μm thick to satisfy the same requirement.
An alternative way to define critical substratum thickness is by comparing not absolute strain values within the ECM, but relative deflections caused by FA traction: at what depth is the in-plane deflection a certain percentage of the in-plane deflection at the FA-ECM interface at the substratum free surface? This can be quantified in terms of FA radius a. By examining Equation 6 (see Methods), we find, for example, that the deflections at depths 0.54a and 3.3a are 50% and 10% of the surface deflection, respectively. We might therefore assign the critical thickness to be 3.3a, or 1.9 μm if we assume a circular FA with a cross-sectional area of 1 μm 2 , since at There is no statistically significant difference in average rupture force of the complex between these two solution pH values. this depth the in-plane deflection is reduced by 90% as compared to the surface deflection. Note that this interpretation of critical thickness, by basis of its definition and assumption of substratum isotropy and constitutive laws, is dependent on FA size, but independent of substratum elastic modulus. An advantage of this second approach is that instead of assigning an arbitrary strain cutoff to define critical thickness, we are comparing at-depth displacements to the displacements exerted by cells at the cell-substratum interface. Continued exploration of these two models of critical thickness and their evaluation by experiment is one focus of current work in our group (in preparation).

Discussion and Conclusion
In order to understand complex processes such as cell adhesion and migration, we must consider chemomechanical features of the materials with which the cell interacts, whether these are synthetic substrata or natural ECM. Computational methods provide an efficient way to test how properties of the extracellular environment can affect time and length scales of focal adhesion formation, stability and force transmission. Simulation techniques such as steered molecular dynamics allow us to gain an atomistic level, mechanistic understanding of how environmental properties such as pH and ECM stiffness can influence the cell-matrix contact and downstream cell properties. With this knowledge, we have the potential to find ways to promote or inhibit certain cell behaviors. Continuum methods such as FEM simulation and analytical mechanics allow us to examine ECM deformation under a variety of conditions, and to understand how cells probe the extracellular environment.
It is important to note that with current computational resources, both of these simulation methods require simplified representations of the cell-matrix interface. In reality, the ECM is comprised of a variety of molecules and also includes embedded growth factors and cytokines. 83 These components would all contribute to microniche properties such as extracellular pH, ECM stiffness and ECM thickness. Furthermore, cell surface integrins can exist in low and high affinity conformations, and associate with each other, intracellular adapter proteins, and other transmembrane receptors during focal adhesion formation. 84 All of these factors can affect cell signaling pathways, which has important consequences for cell behavior. Our simplified models do not incorporate these details, as we neglect all molecular lengthscales in the For t = 0.5 μm, the strain does not decay within the thickness of the ECM, indicating that stress and strain fields will likely be distorted by an underlying rigid material. Note that data if offset for clarity; principal strains for substrata of 3 and 6 μm thickness were commensurate over the distance from the interface u z ; (D) Critical thickness, defined as ECM depth at which FA-induced principal strain is <0.1%, is a function of ECM stiffness, as shown for a total substratum thickness t = 3 μm.
comtinuum-level FEM approach and neglect additional molecular species in the atomistic-level SMD of integrins, but they allow us to systematically study the effects of several important ECM properties without the complication of coupling between different variables. As we gain understanding from such simple models, we as a community can build more complex representations of the cell-matrix interface by incorporating additional details about ECM and FA component proteins as well as cell signaling events.
Our molecular level results for the biotin-streptavidin complex make clear that the forces required to unbind ligand-receptor pairs can show a great deal of variability, over and above the established loading rate dependence. It is important to characterize and understand this variation, because this will determine the overall effect of the multiple ligand-receptor pairs involved in cell adhesion and migration. We demonstrated that RMSD-based clustering analysis is a useful tool for SMD input configuration selection, allowing efficient characterization of this variation. We applied this clustering method to our study of the effect of acidic extracellular pH on integrin α v β 3 -RGD unbinding. We found that an acidic pH of ~5.5 does not significantly change the average rupture force of the complex, and we are currently analyzing possible effects on unbinding pathways and mechanisms.
The biotin-streptavidin simulations also suggest that ECM stiffness can affect the bond lifetime of an ECM ligand and cell surface receptor. If this finding holds true for integrin-ligand systems, this suggests that a complex will exhibit a longer lifetime when the ligand is attached to a more compliant ECM. Here, it is important to note that the overall ECM stiffness in a tissue is more complex than the Hookean spring used in our simulations. To truly represent ECM stiffness, one should take into account the stiffness of individual components as well as the protein network architecture. This is an instance in which models at different length scales could be integrated. For example, a continuum model of the ECM could include subroutine input from a coarse-grained model of viscoelastic/ poroelastic deformation and an atomistic representation of the ligand-receptor complex.
Our continuum level results concerning the lengthscales of FA-induced force transmission have valuable applications toward the design of synthetic substrata for cell biology experiments and biotechnological platforms. Common cell culture materials such as tissue culture plastic are much stiffer than the physiological environments in which cells exist in vivo. It is well established that cell morphology and certain functions are modulated by the mechanical cues of the microenvironment, and in an experiment it is important to properly quantify and perhaps recapitulate the mechanical properties of the physiological environment. Our data show that when synthesizing synthetic substrata of variable stiffness, there is a constraint on the thickness. If the substratum is too thin, the cell traction will induce displacements exceeding the thickness of the substratum such that the cell will "feel" the underlying rigid support. For example, if one intends to culture stem cells on a substratum with stiffness comparable to that of muscle tissue, it is crucial to ensure that this is the mechanical environment actually presented by cell-induced traction over a synthetic surface. Otherwise, the cell response cannot be unambiguously correlated with substrata stiffness. 9 Taken together, our results illustrate how ECM chemomechanics can affect the early and late stage events involved in cell adhesion.
First, ECM stiffness could influence the lifetime of integrin-ligand interactions, which are the physical connections between a cell and the adjacent ECM. This will influence the dynamics of FC/FA formation and stability; for example, on a stiff ECM, the integrinligand binding lifetime could be shorter, resulting in less time for intracellular adhesion proteins to localize. Therefore, one could hypothesize that longer times would be required for the cell to form stable focal adhesions. After stable FAs are formed, cytoskeletal contraction will cause the cell to deform its ECM, and the extent of this deformation will depend on the mechanical properties of the ECM. Cell-induced deformation of ECM proteins such as fibronectin could expose hidden binding sites and induce assembly of fibronectin into its fibrillar form, 36,85 thereby regulating ECM structure and network connectivity. It is also likely that the ECM and/or the cytoskeleton could display nonlinear stress-strain behavior upon large strains. With a nonlinear relationship between stress (σ) and strain (ε), piecewise approximation of dσ/dε at arbitrary Δε is sometimes interpreted as "strain stiffening" or "strain hardening", as it requires an increasingly large change in σ to produce significant increase in ε for strains greater than the threshold value ε p (Fig. 7). If such nonlinearity over physiologically relevant strain magnitudes is significant, then a certain level of cytoskeletal stress could suddenly result in a different level of ECM deformation. This could have interesting consequences for localization of proteins within the ECM, underlying the intricate coupling between mechanics and chemistry in such systems. While a growing body of work has focused chiefly on mechanics within or adjacent to cells, it is imperative to consider the effects of chemical coupling and feedback on the mechanical niche and the cell responses. As these states are inherently coupled in a physiological context, understanding this chemomechanical interplay is an exciting challenge in the ongoing experiments, modeling and simulation of cell adhesion and migration. Figure 7. Schematic illustrating nonlinear stress-strain behavior. There is a linear regime at small strains, but the stress-strain relationship becomes nonlinear at larger strains. For ε > ε p , piecewise approximation of dσ/dε is sometimes interpreted as "strain stiffening" or "strain hardening".

Molecular dynamics simulations of biotin-streptavidin.
Simulation setup and equilibration. The biotin-streptavidin tetramer (PDB ID 1STP 50 ) was simulated using the GROMACS molecular dynamics package, version 3.3 60,61 as described previously. 38 Briefly, the protein was solvated with ions added to provide charge neutrality and to mimic physiological conditions. After steepest descents minimization, unconstrained molecular dynamics simulation over 100 ns was performed to equilibrate the system. Using the protocol developed in ref. 38, we determined that the complex had equilibrated within 15 ns. The equilibrated trajectory (t > 15 ns) was used for further analysis.
Clustering analysis. To efficiently and objectively choose initial configurations for SMD simulations, we used a clustering algorithm to divide the equilibrated portion of the molecular dynamics trajectory into configurational subsets. The GROMACS tool g_cluster was used for single linkage hierarchical clustering. The similarity measure for clustering was the RMSD of residues after least-squares fitting between pairs of configurations. RMSD was calculated with respect to the following groups: (1) full streptavidin tetramer with four bound biotins (RMSD cutoff = 0.75 nm); (2) residues in one streptavidin binding pocket with one bound biotin molecule (RMSD cutoff = 0.35 nm); (3) streptavidin binding pocket residues only (RMSD cutoff = 0.25 nm). The binding pocket was defined as all streptavidin residues with 3.5 Å of biotin in the X-ray diffraction structure (Fig. 8A). 22,45,50 Clusters with at least ten members were chosen for further analysis. The cluster member with the lowest average RMSD with respect to all other cluster members was defined as the centroid of the cluster. Steered molecular dynamics. SMD simulations were performed according to the protocol described in ref. 24. Briefly, one subunit of the streptavidin tetramer was subjected to loading forces through a force transducer (i.e., a Hookean spring). The center of mass of one streptavidin monomer was fixed, but the system was allowed to rotate about the center of mass. The tensile loading direction was defined as the normalized vector between the initial center of mass of the fixed streptavidin subunit and the O2 atom of the biotin bound primarily to that subunit (Fig. 8A). All SMD simulations used the same seed for random initialization of atomic velocities. Transducer spring constants k c ranged from 500 to 5000 kJ mol -1 nm -2 (0. 83-8.3 N/m), while velocities v ⇀ ranged from 0.4 to 10 m/s. Effective loading rates F ' ranged from 0.4 to 11 N/s. The resulting trajectories were analyzed to extract the force exerted by the spring and the reaction coordinate of the ligand as functions of simulation time. The reaction coordinate χ is defined as the distance of the biotin O2 atom from its initial position, χ = [(xx 0 ) 2 + (y -y 0 ) 2 + (z -z 0 ) 2 ] 1/2 . Forces were examined every 200 fs, and rupture force F R was taken as the maximum force recorded during a particular trajectory, which was ultimately followed by dissociation of the complex.
Binding lifetimes were calculated with a least-squares linear regression of F R to F ' to obtain the slope (m) and y-intercept (b) through application of the following adaptation of Bell's model: (1) where k B is Boltzmann's constant, T is absolute temperature, x b is the distance between the bound state and the energetic maximum, 9 Cell F ' = kv ⇀ is the loading rate, and τ is the equilibrium binding lifetime (inverse of the kinetic dissociation rate at equilibrium). From the relationship between F R and F ' given in Equation 1, it follows that Molecular dynamics simulations of integrin-RGD. Simulation setup and equilibration. The extracellular portion of integrin α v β 3 in complex with an RGD ligand (PDB ID 1L5G 62 ) was simulated using GROMACS version 3.3. 60,61 Mn 2+ ions were replaced with Mg 2+ , and only the α propeller, β hybrid and βA domains were simulated for computational efficiency. Particle Mesh Ewald electrostatics were used with a short range interaction cutoff of 0.9 nm. The PRODRG server 63 was used to generate the RGD ligand topology. Multi-Conformation Continuum Electrostatics 64,65 was used to predict accurate pKa values for all the titratable amino acid residues in the system. To simulate an effective acidic pH, all residues with pKa > 6.0 were protonated. To simulate the normal extracellular pH of 7.4, all residues with pKa > 7.9 were protonated. After protonation, the protein was solvated in a box of dimensions 10.235 nm x 11.513 nm x 8.298 nm. For the pH 7.4 system, 28427 water molecules, 91 Na + ions and 81 Clions were added. For the pH 5.5 system, 28428 water molecules, 87 Na + ions and 81 Clions were added. For each pH, a two step steepest descents minimization of the X-ray diffraction structure was performed. In the first step, the integrin, RGD and Mg 2+ ions were held fixed and the maximum force in the system was reduced to less than 2000 kJ mol -1 nm -1 . In the second step, the full system was free to move and the maximum force in the system was reduced to less than 1500 kJ mol -1 nm -1 . After minimization, a 10 ps molecular dynamics simulation was performed with position restraints on the side chains of ARG RGD and ASP 218 . Five simulations were performed for each pH, each with a different seed for random initialization of atomic velocities. In this way, five different input configurations were created for each pH. Next, an 8 ns molecular dynamics simulation was performed with each input configuration to equilibrate the system. During molecular dynamics simulations, center of mass rotation and translation of the receptor were restrained. Systems were considered to be equilibrated when the RMSD of the receptor, ligand, and ions with respect to the initial configuration reached a plateau (usually requiring 2-3 ns).
Clustering analysis. The equilibrated portion of each trajectory was used as input for clustering analysis. RMSD of binding cleft residues after least-squares fitting was used as the similarity measure for clustering. The binding cleft included the 3 Mg 2+ ions in the β subunit as well as all integrin residues that have been reported to interact with the RGD ligand or with the ions 62,66 (Fig. 8B). Clusters with 10+ members were chosen for further analysis. For each trajectory, the RMSD cutoff was chosen in order to have 5-8 clusters of this size. Cluster centroids with ASP RGD -Mg 2+ coordination states similar to those observed in references 32 and 62 were selected as input configurations for SMD simulations.
Steered molecular dynamics. The initial position of the force transducer (spring) in all SMD simulations coincided with the center of mass of the RGD ligand. The centers of mass of the α and β subunits were held fixed. The tensile loading direction was defined as the vector between the initial center of mass of these two fixed points and the center of mass of the RGD ligand. The transducer was defined by stiffness k c = 2492 kJ mol -1 nm -2 (4.14 N/m) with v ⇀ = 2 m/s (Fig. 8B); this corresponds to F ' = 8.28 N/s. All SMD simulations used the same seed for random initialization of atomic velocities. Analysis was performed with the procedure described above for biotin-streptavidin.
Continuum simulations of FA traction against the ECM. A three-dimensional finite element model of a FA-ECM complex was constructed using 8-noded brick elements. Due to symmetry, only one quarter of the system was modeled to obtain the stress field evolution within the FA and the underlying substratum during traction application. After performing mesh sensitivity analysis, the optimized model consisted of 16,434 nodes and 14,355 elements. The mesh generator TRUEGRID (XYZ Scientific, Livermore, CA) was used to mesh the FA-ECM complex, and the finite element analysis was performed using the general purpose non-linear finite element analysis code ABAQUS (HBK, Providence, RI). Figure 9 shows the FE mesh with substratum thickness t = 6 μm. The FA was modeled with a square geometry, as initial studies indicated negligible effects on substratum deformation for more complex circular or ovoid FAs. FA size and thickness were based on estimates from immunofluoresent staining of FAs in human fibroblasts. 67 The FA  was idealized as linear elastic, with E = 10 kPa and Poisson's ratio ν = 0. 3,8,68 while the ECM was modeled as a nearly incompressible linear elastic solid (ν = 0.49) with E = 18, 180 and 1800 kPa, each with ECM thickness t = 6 μm, 3 μm and 0.5 μm. The lowest E (18 kPa) is representative of the elastomeric substrata used by Balaban et al., in experiments designed to measure FA traction magnitude, 67 and the highest E (1800 kPa) is in the range of PEM substrata stiffness that have been employed by our group and others in the study of mechanosensitive adhesion of tissue cells. 14 The consideration of three elastic moduli and three thicknesses allows decoupling of material and geometric parameters, as well as a comparison of traction generated when the effective stiffness of the focal adhesion and the adjacent ECM are well-matched. Shear stress τ was applied to the intracellular surface of the FA complex with τ = 5 kPa (the experimentally observed traction exerted at FAs by fibroblasts in ref. 67), and the resulting von Mises effective stresses generated through the compliant FA were computed as a function of depth into the adjacent ECM. Note that von Mises stress is defined as a coordinate-system independent combination of principal stresses σ eff that results in permanent deformation of the material when σ eff exceeds the uniaxial yield stress of the material σ y : where σ i = principal (normal) stress. Analytical model of FA-induced deformation of the ECM. We have also analytically modeled substratum deformation caused by FA traction. We assume the substratum to be an elastically isotropic and homogeneous half space, and to deform incompressibly (Poisson's ratio ν = 0.5) within the linear elastic regime. We idealize the FA as a circular area S of radius a over which a uniform shear stress τ is exerted. The horizontal displacement in the direction of τ at a depth z under the center of the FA is (4) where we are integrating the Green's tensor G for a tangential point load on the surface of a half space. This Green's tensor for an incompressible material 69 is (5) where , E is the Young's elastic modulus of the substratum, and the x-axis corresponds to the direction of FA-applied shear stress τ. By integrating over the circular region S, we obtain the displacement function (6) This result is valid for substrata that are sufficiently stiff that the loading geometry is not altered by the resulting deformation. We can estimate a lower limit of suitable substratum elastic modulus by examining the change in deflection with depth: (7) This value represents the maximum slope of an originally vertical, imaginary line located at the center of the FA. To maintain the tangential component of the shear stress within 20% of its original value, we require the angle of the line to be less than 0.64 radians with respect to vertical, which can be shown through trigonometry to be equivalent to requiring that E > 4τ. Assuming a FA shear stress τ = 5 kPa, Equation 6 is therefore reasonably accurate for a linear elastic substratum stiffness E ≥ 20 kPa.