Cells interacting through an extracellular matrix (ECM) exhibit emergent behaviors resulting from collective intercellular interaction. In wound healing and tissue development, characteristic compaction of ECM gel is induced by multiple cells that generate tensions in the ECM fibers and coordinate their actions with other cells. Computational prediction of collective cell-ECM interaction based on first principles is highly complex especially as the number of cells increase. Here, we introduce a computationally-efficient method for predicting nonlinear behaviors of multiple cells interacting mechanically through a 3-D ECM fiber network. The key enabling technique is superposition of single cell computational models to predict multicellular behaviors. While cell-ECM interactions are highly nonlinear, they can be linearized accurately with a unique method, termed Dual-Faceted Linearization. This method recasts the original nonlinear dynamics in an augmented space where the system behaves more linearly. The independent state variables are augmented by combining auxiliary variables that inform nonlinear elements involved in the system. This computational method involves a) expressing the original nonlinear state equations with two sets of linear dynamic equations b) reducing the order of the augmented linear system via principal component analysis and c) superposing individual single cell-ECM dynamics to predict collective behaviors of multiple cells. The method is computationally efficient compared to original nonlinear dynamic simulation and accurate compared to traditional Taylor expansion linearization. Furthermore, we reproduce reported experimental results of multi-cell induced ECM compaction.
Cell-induced compaction of fibrous extracellular matrix (ECM) is an important mechanism for numerous processes such as wound healing and tissue development [1–3]. During wound healing, for example, traction forces exerted by fibroblasts and myofibroblasts result in ECM compaction at the site of injury [2,3]. In vitro experiments using cell-populated collagen gel reveal global compaction of the matrix as a result of cooperative effect of multiple cells at the boundaries as well as propagation through the bulk [4–6]. Furthermore, matrix densification is observed in the regions around [7] and in-between cells. Here we examine the mechanical aspect of intercellular communication through the ECM and how contractile cells can induce emergent mechanical changes leading to matrix compaction. From a simplified mechanics point of view, compaction results when the traction forces exerted by the contractile cells embedded within the ECM overcome the resistive forces of the ECM structure, including viscoelastic forces and elastic energy forces. As a result the matrix is deformed from its original stress-free state and the elastic modulus increases [4–7]. In reality, the compaction process is far more complex. The ECM forms a network of cross-linked fibers that is highly nonlinear and intricate, but is critical for predicting large compaction and long-range transmission of forces [4]. As a large deformation is induced by contractile cells, the standard linear mechanics model yields substantial errors. The ECM fiber network is anisotropic and causes irreversible deformations as a large compaction takes place. This prominent nonlinearity prohibits use of simple methods for predicting the ECM compaction by a multitude of cells. In addition, cells can internally modulate their state in response to local mechanical stresses within the ECM, which influences cell polarity, contractility, stiffness and strength of focal adhesion’s [8–11]. These cell properties are highly nonlinear and complex. Consideration of these nonlinear physical and physiological properties involved in the cell-ECM mechanics often result in differential equations that are intractably complex due to high-dimensional, nonlinear coupled dynamics. Many in silico modeling approaches in the areas of wound healing and fibrotic disease have helped elucidate and explore the underlying phenomena involved in cell-induced ECM compaction, and have been used to supplement in vitro experiments for fast and inexpensive methods of evaluation. Approaches in previous works include: i) a hybrid continuum-discrete framework consisting of the macroscopic finite element domain and local microscopic fiber network [12], ii) rule-based models with deformable cells and ECM fibers to explain matrix remodeling and durotaxis [13,14], iii) a discrete fiber model of cell populated fibrous matrix [15], and iv) continuum models of ECM gel compaction [7,16,17]. Even though these works provide many insights, they also simplify the ECM gel compaction mechanism by: a) 2-D representation of a 3-D system, b) exclusion of intracellular mechanics such as mechanobiology of actin stress fibers, focal adhesions, and remodeling of cellular and nuclear membranes, and c) consideration of linear elastic spring model of ECM fibers without including the viscoelastic nature of the fibers. Consequently, these prior models abstract detailed cell-ECM interactions, resulting in limitations to understanding how these interactions enable characteristic gel compaction. In addition, models examining the complex dynamics surrounding cell morphology, contractility, and polarity based on finite element methods do exist for 2 dimensional cases [10]. And, focal adhesion-stress fiber dynamics have been modeled for 2-D PDMS substrates based on non-equilibrium thermodynamics [11]. In the current work, the ECM is modeled as a 3-D cross-linked network of discrete, viscoelastic fibers, and detailed mechanistic cell dynamics, including focal adhesion dynamics, cytoskeleton remodeling, actin motor activity and lamellipodia protrusion, are derived from basic principles. The resultant model is computationally complex, especially for a larger number of cells. The governing differential equations are highly nonlinear, coupled, and of high dimension. Here, we solved this difficulty by introducing a methodology having its disciplinary basis spanned in system dynamics, machine learning, and statistics. It is known that a nonlinear system can behave more linearly when recast in a larger space [18]. In our approach, the original nonlinear dynamics derived from physical and physiological principles are recast in an enlarged state space by augmenting independent state variables with auxiliary variables that inform input-output characteristics of the nonlinear elements involved in the system. Once recast in the augmented space, the nonlinear system can be represented as an augmented set of linear dynamic equations. The linear representation facilitates model reduction using latent variable analysis, which can be shown is difficult to apply to highly nonlinear systems [19–22]. Furthermore, linearization in the augmented space allows for superposition of multiple subsystems. In the current work, collective behaviors of multiple cells are predicted via superposition of single cell subsystems through the linearization in the augmented state space. The proposed methodology is general, and is applicable to a broader class of problems where large-scale, collective behaviors must be predicted while retaining sufficient mechanistic details.
We construct a computational model for predicting cell-mediated gel compaction by multiple (ncell) cells having a uniform phenotype and interacting through a surrounding 3-D ECM fiber network. The ECM is modeled as a network of many fibers connected at a large number of nodes (Ne ≈ 2000), whereas each cell is represented with a mesh structure consisting of multiple nodes (Nc ≈ 200) which forms the cell outer membrane (see Fig 1A). The cell outer membrane deforms and gains traction as the nodes on the membrane bond to the nodes of the surrounding ECM fiber network and form focal adhesions, which occur when bonding molecules (or integrins) on the cell membrane bind to ligands on ECM. Consider the i-th outer membrane node of the k-th cell with three dimensional spatial coordinates x i c, k ∈ ℜ 3 × 1 (See Fig 1B). The forces acting on it include the cell’s cortical tension force and elastic energy force (collectively denoted as F C o r t - E l a s, i c, k ∈ ℜ 3 × 1), focal adhesion force (denoted as F F A, i c, k ∈ ℜ 3 × 1), lamellipodium force (F L, i c, k ∈ ℜ 3 × 1), and frictional damping force (F D a m p, i c, k ∈ ℜ 3 × 1) [23,24]. Assuming that the mass of the node is negligibly small and the damping force is given by F D a m p, i c, k = - D c d x i c, k / d t, where Dc is damping constant, the equation of motion is given by: F C o r t - E l a s, i c, k + F F A. i c, k + F L, i c, k - D c d x i c, k d t = 0 i = 1, ⋯, N c, k = 1, ⋯, n c e l l (1) The generation of lamellipodium force pertains to the polarity of the cell. Namely, lamellipodia extend in a particular direction of the cell determined by the cell’s polarity [23–26]. The cell polarity and the lamellipodium forces can be treated as a cell’s decision or, in the system dynamics terminology, control inputs. Let x c, k = (x 1 c, k T ⋯ x N c c, k T) T ∈ ℜ 3 N c × 1 be a vector containing the 3-D coordinates of all the cell membrane nodes. Here the superscript in XT represents the transpose of matrix or vector X. The above equation of motion can be written collectively as: d x c, k d t = W C E c F C o r t - E l a s c, k + W F A c F F A c, k + L c u k k = 1, ⋯, n c e l l (2) where F C o r t - E l a s c, k ∈ ℜ 3 N C × 1 is a vector comprising cortical tension and elastic energy forces for all the cell nodes (i = 1, ⋯, NC), F F A c, k ∈ ℜ 3 N C × 1 is a vector of focal adhesion forces at all the cell nodes, u k ∈ ℜ 3 N C × 1 is an input vector containing all the lamellipodium forces (F L, i c, k), and W C E c, W F A c and Lc are constant matrices of consistent dimensions. The equation of motion of the surrounding ECM fiber network can be represented in a similar manner. The forces acting on the j-th node of the fiber network are the elastic energy forces, including both lateral restoring forces and the one associated with bending moments, (F E l a s, j e ∈ ℜ 3 × 1), focal adhesion forces from the shared attachment with the cell (F F A, j e ∈ ℜ 3 × 1) and damping forces (F D a m p, j e ∈ ℜ 3 × 1) [23–26]. The equation of motion can be written as: F E l a s, j e + F F A, j e - D e d x j e d t = 0, j = 1, ⋯, N e (3) Let x e = (x 1 e T ⋯ x N e e T) T ∈ ℜ 3 N e × 1 be a vector containing the 3-D coordinates of all the ECM nodes. Then Eq 3 can be written as: d x e d t = W E l a s e F E l a s e + W F A e F F A e (4) The ECM elastic energy force is a nonlinear function of ECM coordinates xe. The cortical tension and elastic energy force of the k-th cell is a nonlinear function of its membrane coordinates xc, k. Here xe and xc, k are independent state variables of the multi-cell ECM system. F E l a s e = F E l a s e (x e), F C o r t - E l a s c, k = F C o r t - E l a s c, k (x c, k) k = 1, ⋯, n c e l l (5) The focal adhesion force is modeled as a stochastic binding process between nodes on the cell membrane and those on the ECM. Using Monte Carlo simulations it has been found that focal adhesion forces can be approximated to a nonlinear algebraic function of cell membrane and ECM nodes as well as the biochemical parameters involved in integrin-ligand binding [23,24]. F F A c, k = F F A c, k (x c, k, x e), k = 1, ⋯, n c e l l (6) Assuming that no two cells bind to the same ECM node, we can find that the focal adhesion force of the i-th membrane node of the k-th cell attached to the j-th ECM node must satisfy: F F A, i c, k + F F A, j e = 0 (7) Namely, F F A, i c, k and F F A, j e have the same magnitude with the opposite signs. Therefore, all the focal adhesion forces of the k-th cell can be mapped to the corresponding ECM nodes. Collectively, the focal adhesion forces of all the nodes within the ECM may be written as: F F A e = ∑ k = 1 n c e l l P m a p k F F A c, k (8) where P m a p k ∈ ℜ 3 N e × 3 N c is a parameter matrix (consisting of either 0 or -1 elements) that maps the membrane focal adhesion forces of the individual cells (F F A c, k, k = 1, …, n c e l l) to the corresponding ECM focal adhesion forces (F F A e). The focal adhesion connections between the membrane nodes and ECM nodes change over time as the cell membrane deforms, gains traction and generates lamellipodial protrusions. Therefore, the mapping matrix P m a p k is updated at each time step. Details on the formation and structure of P m a p k are given in the Methods Section. The ncell cells interact with each other through the surrounding ECM by generating focal adhesion forces, which propagate through the ECM fiber network and influence the other cells. The resultant collective behavior of the multiple cells is complex due to coupled, nonlinear dynamics. Although the governing equations derived above are rigorous and based on basic principles, they are complex and can become computationally expensive as the number of cells increases. Computational complexity is a key challenge in predicting collective behaviors of multiple cells. The number of state variables for the given system is 3Ne + 3Nc ncell, which is on the order of 7,000 for ncell = 2 and 9000 for ncell = 5. We aim to transform the governing equations into a linear latent variable representation in order to considerably reduce the number of state variables but also facilitate prediction of collective behaviors of the multiple cells through superposition of individual cell dynamics. Model reduction is a challenging problem particularly for highly nonlinear, dynamical systems [21,22,27–29], as in the presented problem of collective behaviors of multiple cells within an ECM. If the system is linear or near linear, model reduction is more amenable and simple methods, such as Principal Component Analysis and Partial Least Squares, can reduce dimensionality. Here, we propose a unique linearization method, termed Dual-Faceted (DF) Linearization, and then apply a model reduction method to the linearized model. In DF Linearization, we represent the nonlinear dynamical system in an augmented space consisting of independent state variables (xe and xc, k) and nonlinear forces (F E l a s e, F C o r t - E l a s c, k, F F A c, k) as the additional variables, termed auxiliary variables. Standard linearization, such as Taylor series expansion, is limited in accuracy, which may be valid only in the vicinity of a reference point. In DF Linearization, instead of taking “algebraic” linearization of these nonlinear terms, we consider “dynamic” linearization by representing their dynamic transitions using linear regressions. Before formally introducing DF Linearization, let us consider a simple example that delineates the basic principle of the method. Suppose that the system consists of one spring and a damping element with negligibly small mass, F - D x ˙ = 0. If the spring is a linear spring, F = kx, there is absolutely no difference between the equation in terms of state variable x, x ˙ = (k / D) x, and the one in force F, F ˙ = (k / D) F. However, it is not the case if the spring is nonlinear, for example a hard spring: F = ax + bx3 where a > 0, b > 0. Representing the differential equation in two variables, one with the state variable x and the other with the auxiliary variable F, provide different equations. d x d t = a D x + b D x 3 (9) d F d t = 1 D (a + 3 b (g (F) ) 2) F (10) where x = g (F) s the inverse function of F = ax + bx3. Both equations represent the same nonlinear system, yet the expressions are different, hence Dual-Faceted representations. Linearizing these differential equations lead to two linear differential equations viewed from the augmented space. Note that Eq 9 can be represented as a linear equation by using both state and auxiliary variables: d x d t = W F F (11) where WF = 1/D. The augmented state Eq 10 can be approximated to a linear regression: d F d t ≃ S x x + S F F (12) where Sx, SF are regression parameters. The expression given by Eq 12 differs from the one based on the first order Taylor expansion (or “algebraic” linearization) which yields: F (x) ≃ F (x ¯) + d F d x Δ x (13) Furthermore, if we evaluate the derivative J (x) ≡dF/dx at a particular point, J ¯ (x ¯), and then use Eq 13 to express the augmented state equation, it reduces to F ˙ = J ¯ x ˙. This implies that F ˙ and x ˙ are proportional. Using an “algebraic” linearization yields a differential equation representing the transition of F that is collinear to the one representing the transition of x, and thereby an auxiliary state equation would not provide any new information. Conversely, the regression model in Eq 12 provides us with a diverse view of the original nonlinear system, thus providing a richer representation of the nonlinearity than the standard first order Taylor expansion. Applying the above principle of Dual Faceted Linearization to our problem, we note that the original state equations governing the transition of the independent state variables, 2 and 4, are linear in the augmented state space. All we need is to obtain the transition of the auxiliary variables. Let the regression of the dynamic transition of auxiliary variable F E l a s e, be expressed as: d F E l a s e d t ≈ R x e x e + R F E l a s e F E l a s e (14) where R x e, R F E l a s e ∈ ℜ 3 N e × 3 N e are parameter matrices. If an “algebraic” linearization using the Jacobian J ¯ = ∂ F E l a s e / ∂ x e | x ¯ e was utilized, the above equation would be: d F E l a s e / d t = J ¯ · d x e / d t. This state transition equation through “algebraic” linearization is equivalent to one of the original independent state Eq 4 because d F E l a s e / d t and d x e / d t are collinear within this formulation which renders it redundant. In contrast, the state transition equation presented in Eq 14 is not collinear, providing a diverse facet of the nonlinear system. Similarly, for the auxiliary variables F C o r t - E l a s c, k, F F A c, k, let the regression equations be written as: d F C o r t - E l a s c, k d t ≃ Q x c x c, k + Q F C E c F C o r t - E l a s c, k + Q u u k (15) d F F A c, k d t ≃ H x c x c, k + H x e x e + H F F A c F F A c, k + H u u k (16) where Q ⋆ ⋆, H ⋆ ⋆ (⋆ -corresponding subscript and superscript) are parameter matrices with consistent dimensions. The high-dimensional parameter matrices (R ⋆ ⋆, Q ⋆ ⋆, H ⋆ ⋆) do not need to be determined explicitly as discussed in the subsequent sections. DF Linearization represents a nonlinear dynamical system with two sets of differential equations. One set is the original state equations governing the transition of the independent state variables and the other set is the regression of the dynamics of auxiliary variables. The original state Eqs 2 and 4, are apparently linear in terms of the auxiliary variables and inputs. In these equations, all the forces acting on each node sum to zero. These are linear expressions when the nonlinear forces are treated as auxiliary variables. In addition, the auxiliary state transitions (Eqs 14–16) are given by linear regressions in the augmented space. Therefore, both differential equations are linear. The two linear differential equations represent different (or dual) facets of the original nonlinear system viewed from the augmented space, thus providing a richer representation of the nonlinearity. Now that the original nonlinear system has been represented as a linear dynamical system in the augmented space, we can apply a latent variable modeling method to reduce model order. Represented in the augmented space, the differential equations may contain similar modes, or some variables are close to collinear. These similar modes and collinear variables can be eliminated by model order reduction methods. Let ζ c, k be the augmented variable vector containing membrane node coordinates and forces of the k-th cell. ζ c, k = (x c, k F C o r t - E l a s c, k F F A c, k) ∈ ℜ 9 N c × 1 (17) Here uk (the cell’s lamellipodial force) is treated as input variables that are excluded from the augmented state space. Similarly, let ζe be the augmented variable vector containing ECM node coordinates and forces: ζ e = (x e F E l a s e) ∈ ℜ 6 N e × 1 (18) Focal adhesion forces F F A e are determined by the individual cells by Eq 8 and, thereby, excluded from the augmented space of the ECM. We apply latent space analysis to vectors ζc, k and ζe, respectively. First we generate data by using Eqs 2 and 4–7. Computation of the nonlinear state equations is amenable for a single cell interacting with ECM. A data set can be created by simulating those nonlinear equations by placing a single cell at diverse locations, i. e. repeating the simulation with different initial conditions. Let C ζ ζ c, C ζ ζ e be the covariance matrices of simulation data sets of augmented states ζc, k and ζe, respectively. Each covariance matrix contains both independent state and auxiliary variables, where the latter is nonlinear functions of the former. If auxuliary variables were linear functions of the state variables, then the rank of the covariance matrix would be equal to the number of independent state variables. However due to the nonlinearity, the rank is higher. Details on the formation of C ζ ζ c, C ζ ζ e are given in the Methods Section. The covariance analysis also reveals that the system represented in the augmented space contains many components that may be negligibly small. Using Principal Component Analysis, the original data of ζc, k and ζe can be represented with latent variables of truncated dimension mc ≪ 9Nc and me ≪ 6Ne, respectively [29]: ζ c, k = (V x c V F C E c V F F A c) ︸ V c z c, k, ζ e = (V x e V F E l a s e) ︸ V e z e (19) where matrices V c ∈ ℜ 9 N c × m and V e ∈ ℜ 6 N e × m are orthogonal matrices comprising the eigenvectors of the covariance matrices, and z c, k = (V x c V F C E c V F F A c) T (x c, k F C o r t - E l a s c, k F F A c, k) ∈ ℜ m c × 1 z e = (V x e V F E l a s e) T (x e F E l a s e) ∈ ℜ m e × 1 (20) Differentiating the latent variable state vector zc, k and substituting Eqs 2,15,16 and 19 yields: d z c, k d t = V x c T d x c, k d t + V F C E c T d F C o r t - E l a s c, k d t + V F F A c T d F F A c, k d t = A z c, k + B u k + C z e, k = 1, ⋯, n c e l l (21) where: A = V x c T (W C E c V F C E c + W F A c V F F A c) + V F C E c T (Q x c V x c + Q F C E c V F C E c) + V F F A c T (H x c V x c + H F F A c V F F A c) B = V x c T L c + V F C E c T Q u + V F F A c T H u C = V F F A c T H x e V x e (22) Eq 21 provides the latent variable state equation of the k-th cell interacting with the ECM. Given the latent variable state of ECM ze and the input uk reflecting the cell’s decision, the transition of the cell’s latent variable state is determined locally without directly including the states of the other cells. Cells interact indirectly through the strain field created by other cells over the ECM fiber network. The ECM dynamics can be represented in the latent variable space spanned by Ve. Differentiating the latent variable state vector and substituting Eqs 4,14 and 19 yield: d z e d t = V x e T d x e d t + V F E l a s e T d F E l a s e d t = G z e + ∑ k = 1 n c e l l D k z c, k (23) where: G = V x e T W E l a s e V F E l a s e + V F E l a s e T (R x e V x e + R F E l a s e V F E l a s e) D k = V x e T W F A e P m a p k V F F A c (24) Fig 2 shows numerical examples of the DF linearization and subsequent latent variable transformation in reproducing accurate cell morphologies of the original nonlinear computational model over time. Remarkably, the DF linearized latent variable model can correctly reconstruct the complex cell membrane topology with m = mc + me = 50 + 50 = 100 total latent variables. Fig 2C quantifies the root mean square error and computation time as a function of latent variable model dimension. As can be seen, the computation time for the latent variable cell-ECM system increases with increased latent variables while the root mean square error decreases. Conversely, the standard Taylor expansion linearization is not capable of representing cell morphologies without marked error which is quantified in Fig 2B. We also compare our DF linearization approach to a more sophisticated method for approximation of nonlinear systems termed trajectory piece-wise linear (TPWL) [30,31]. This method uses collection of (algebraic) linearizations of the original nonlinear system about suitably selected states to approximate the nonlinear system. As can be seen from Fig 2B, although the TPWL method yields lower error than the first order taylor expansion, our DF linearized model still leads to significantly lower prediction error. This is because in DF linearization, instead of taking “algebraic” linearization of nonlinear terms, we consider “dynamic” linearization by representing their dynamic transitions using linear regressions. These results demonstrate the effectiveness of DF linearization and model reduction in reconstructing simulations from a high dimensional complex nonlinear model. This latent space model provides not only a low-dimensional structure for efficient computation, but also contains natural insights into the interactions among the multiple cells. Fig 3A shows the dynamic interactions in block diagram form based on Eqs 21 and 23. The ECM changes its latent variable state ze with the autoregressive feedback through matrix G as well as with the forward path that collects the latent variable states of all the individual cells, as shown by the summing junction Σ. Each single cell changes its latent variable state zc, k (k = 1, …, ncell) with autoregressive feedback through A as well as with two forward paths. The first path (fed through C) represents global feedback from the ECM (ze). The second path (fed through B) represents the updated lamellipodial forces uk determined from ECM state ze. The lamellipodial forces can be thought of as the individual cell’s decision based on its position and updated ECM properties as explained more in the following section. The actions taken by all the cells are integrated into the global ECM state transition, which is fed back to the individual cells. Therefore, each cell is connected to other cells through the global feedback of the ECM latent variable state ze. Fig 3A manifests the control-theoretic interpretation of multiple cells interacting through ECM. Since the system is represented in a lower dimensional space, the high dimensional regression coefficient matrices (R ⋆ ⋆, Q ⋆ ⋆, H ⋆ ⋆) are not computed explicitly. Instead, the lower dimension coefficient matrices A, B, C, G are computed from numerical simulation data that can be transformed into latent variable space. Details are given in the Methods Section. As discussed previously, input vector uk pertains to the lamellipodial forces generated at each membrane node within the leading edge of the cell. The cell continuously updates its lamellipodial protrusions depending on the orientation of the leading edge as the cell’s polarization (or polarity) changes. The polarity of a cell is important to determine the orientation of the leading edge and is influenced by the direction of local maximum stiffness in the ECM [23–26]. Here we aim to extend the dynamics model of cell polarity developed in [23–25] to predict the formation of lamellipodia. Let d P o l k ∈ ℜ 3 × 1 be a 3-dimensional unit vector indicating the direction of polarity in the k-th cell and d M a x - S t i f f e, k ∈ ℜ 3 × 1 be a 3-dimensional unit vector pointing in the direction of the maximum stiffness of ECM in the vicinity of the k-th cell’s current location. According to [25], the cell polarity rotates dynamically in response to ECM’s local stiffness in such a way that the polarity vector may align with the direction of the maximum stiffness: d d P o l k d t = κ d P o l k × (d M a x - S t i f f e, k × d P o l k) (25) where × indicates vector product, and κ is a scalar parameter. Fig 3B illustrates this relationship. The polarity vector d P o l k tends to align with the maximum stiffness direction, d M a x - S t i f f e, k. The leading edge of the cell is indicated by a right circular cone with apex angle 2 α L k having its centerline aligned with the polarity direction. The membrane nodes of the k-th cell within the cone have nonzero lamellipodial forces (F L, i c, k ≠ 0). Membrane nodes outside this cone have zero lamellipodial forces (F L, i c, k = 0). The direction of maximum ECM stiffness d M a x - S t i f f e, k depends on the stress field within the ECM, which pertains to the latent state vector of ECM ze. The details are given in Appendix C in S1 Text. Compaction of the ECM by the collective efforts of multiple cells is numerically analyzed based on the model reduction and superposition of the nonlinear cell-ECM dynamics via DF Linearization. We first consider the case where two cells placed 30 μm apart are embedded in a 3-D cylindrical ECM that measures 40 μm in diameter and 100 μm in length as seen in Fig 4A. The boundary conditions of the ECM fiber network are set such that the two flat planes on sides are fixed to space (constrained), while the curved surface surrounding the ECM is kept free (unconstrained). The volume of the cylindrical ECM shrinks over time from its initial unstressed state as the cells interact with the surrounding ECM. To quantify the spatiotemporal compaction process, the original ECM cylinder is segmented into 10 slices of 10 μm thickness along its longitudinal axis as shown in Fig 4B. The volumetric changes to the individual slices are plotted in Fig 4C. The prediction of decreased cell volume by the latent variable superposition model (blue) agrees well with the ground-truth, full-scale nonlinear simulation results (green). This is further verified by the corresponding cross-sectional images of the 2-cell cylindrical ECM simulations in Fig 4C. The polarity directions of both cells (shown by red arrows initially pointing in arbitrary directions) shift to point inward, indicating that larger stresses are detected in the area between the cells. A video of the simulation comparison is shown in S1 Video. The proposed model is able to reproduce collective behaviors of multiple cells causing the characteristic compaction of ECM gel, which is not observed for single isolated cell models. This is further verified in Fig 4D which compares the ECM compaction results between single cell and two cell models. As can be seen, the single cell model predicts more localized shrinkage of the ECM volume whereas the two cell model shows more global shrinkage extended to within the region between the cells. A video of the simulation comparison is shown in S2 Video. Fig 4D suggests the presence of more than one cell is necessary for the pronounced ECM compaction leading to emergent changes within the ECM. However, the emergence of pronounced compaction entails not only plurality of cells but proper cell spacing. Fig 5A shows that, as the spacing between cells increases, compaction is less pronounced between them, indicating decreased interaction and integration of cell induced propagated forces. This is summarized in Fig 5B which quantifies the average ECM elastic force in-between cells against cell spacing. From Fig 5B, we see that the average ECM elastic force in-between cells spaced at 100 μm is an order of magnitude less than that of the cells spaced at 30 μm. A video of this simulation is shown in S3 Video. The above computational results shown in Figs 4 and 5 verify that the proposed method can capture collective behaviors of multiple cells. The verification was made by comparing the reduced-order superposition model using DF Linearization against the full-scale, nonlinear model. To further verify the capability of the reduced-order superposition model, a comparison is also made against in vitro experimental data of ECM compaction by a larger number of cells. As shown in Fig 5C and 5D, the computational model successfully reproduces the in vitro experiment conducted by Fernandez, et al [5], in which a heterogeneous planar distribution of MC3T3-E1 osteoblasts were plated in 3-D rectangular prism collagen gel of 50 μm height, 100μm width and 250 μm length. The boundary conditions of ECM in the computational model were set to be consistent with experimental conditions in [5]. The multi-cell latent variable simulation is able to predict the characteristics of the ECM over time. Whereas the group of 5 cells at the left edge exhibit anisotropic contraction of the ECM at the boundary, the single isolated cell at the right edge barely contracts the gel. Fig 5C compares the isolated cell to the group of cells in terms of maximum edge displacement. The isolated cell’s ECM edge displacement is so small that five times its displacement is still substantially lower than the displacement of the group of cells. A video of this simulation is shown in the S4 Video. The presented method for predicting collective behaviors of cell-mediated ECM gel compaction is scalable. Since the individual cell-ECM interactions are local computations, as given by Eq 21, the computational complexity does not increase exponentially, although the number of cells increases.
The nonlinear state Eqs 2 and 4–6 were computed with custom C-code based on references [23,24]. The computation of a single isolated cell embedded in the cylindrical ECM took approximately 24 hours for a single simulation of physical time T ≈ 3600 seconds with sampling interval of 1 second. The simulation was repeated for N ≈ 10 times at different initial cell locations each time. The simulations with a single isolated cell embedded in the large rectangular ECM for reproducing the experimental result were run for approximately 5 days (120 hours). The physical time of simulation was T ≈ 3600. The simulation was repeated for N ≈ 10 for various initial locations of a cell. For each simulation, the total number of sample points for all the variables was over 5,000,000. The number of sample points was 5 × 107. The computation was performed on Intel Xeon CPU E5-2687W @ 3. 10 GHz (2 processors) with 32 logical cores. More details on the formulation of the full-scale nonlinear state equations are summarized in Appendix A in S1 Text. We create the covariance matrix using simulated training data. The training data consists of 3,600 time points of both state and auxiliary variables of a cell embedded in an ECM environment. The simulation is repeated N = 10 times, each time with the cell embedded in distinct locations within the ECM. Then the data covariance matrices may be formed: C ζ ζ c = 1 K · N · T ∑ k = 1 K ∑ n = 1 N ∑ t = 1 T ζ ˜ c, k, n (t) ζ ˜ c, k, n (t) T C ζ ζ e = 1 N · T ∑ n = 1 N ∑ t = 1 T ζ ˜ e, n (t) ζ ˜ e, n (t) T (26) where ζ ˜ c, n, k (t) represents the mean centered t-th time sample (of augmented variable vector ζc, k) for the k-th cell (here K = 1 or 2) embedded within for the ECM at the n-th simulation, and ζ ˜ e, n (t) represents the mean centered t-th time sample of the augmented variable vector ζe in the n-th simulation. By performing eigen-decomposition on the covariance matrices we obtain the orthogonal matrix Vc, Ve comprised the eigenvectors of the data covariance matrix: C ζ ζ c ≈ V c Λ c V c T ∈ ℜ m c × m c C ζ ζ e ≈ V e Λ e V e T ∈ ℜ m e × m e (27) where Λc, Λe are diagonal matrices containing the largest mc and me eigenvalues of the covariance matrices, respectively. It is important to check whether the covariance matrices contain sufficiently rich data, and their first mc and me components are sufficient to capture the cell-ECM dynamics at any cell location within the ECM. Standard techniques can be applied to validate the data and truncation of components [29]. With these, the ECM dynamics of a cell embedded within the ECM at an arbitrary location will be well represented in linear latent variable space which is critical to the success of the method. Covariance matrix calculation were conducted by using Matlab. The following outlines the steps to compute coefficient matrices A, B, C, G by Eqs 21 and 23: Parameter matrices A ∈ ℜ m c × m c, B ∈ ℜ m c × 3 N c, C ∈ ℜ m c × m e, G ∈ ℜ m e × m e are substantially lower in dimension than the regression coefficient matrices R * *, Q * *, H * * given in Eqs 14,15 and 16. Therefore, fewer data points allow us to determine these parameter matrices in the latent variable space. It should be noted that matrices Dk’s are of high dimension, but are not computed with regression since they consist of known matrices as shown in Eq 24. Matlab was used for estimation of parameter matrices and subsequent computations of the latent variable model. 3-D visualization of simulation data was conducted using Tecplot 360. When a focal adhesion is formed between the i-th node of the k-th cell and the j-th node of ECM, the two focal adhesion forces sum to zero, as described previously (F F A, i c, k + F F A, j e = 0 where F F A, i c, k, F F A, j e ∈ ℜ 3 × 1). Representing this relationship in terms of the collective focal adhesion force vectors, F F A c, k ∈ ℜ 3 N c × 1 and F F A e ∈ ℜ 3 N e × 1, requires a matrix P m a p k ∈ ℜ 3 N e × 3 N c. Let F F A e, k be the forces acting on the ECM nodes caused by focal adhesion between the k-th cell and the ECM nodes. This can be written as i ↓ F F A e, k = (⋮ F F A, j e, k ⋮ ⋮) = j → (⋯ ⋯ ⋯ ⋯ ⋯ ⋯ - I 3 × 3 ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯) ︸ P m a p k (⋮ ⋮ F F A, i c, k ⋮) (28) where I3×3 is the 3-dimensional identity matrix. Obtaining this mapping matrix P m a p k for all the cells, the complete focal adhesion forces in the ECM can be expressed in relation to the cell’s focal adhesion forces. F F A e = ∑ k = 1 K F F A e, k = ∑ k = 1 K P m a p k F F A c, k (29) As previously mentioned, the focal adhesion connections between the cell membrane nodes and ECM nodes can vary over time as the cell membrane deforms, gains traction, and generates lamellipodial protrusions. Therefore, P m a p k is updated to reflect the new focal adhesion attachments and detachments at each time step. The original nonlinear computational model has developed a functional relationship between the focal adhesion force, number of integrins, and distance between the membrane and ECM node (see details supplementary materials Appendix A in S1 Text). In the presented framework, the change in focal adhesion attachments can be derived from simulated training of the nonlinear computational model.
The collective ECM compaction by multiple cells is predicted through superposition of individual cells’ contributions in latent variable space. This is made possible by DF Linearization, latent variable transformation and subsequent superposition of single-cell models to predict the collective behavior among multiple cells. As shown in Fig 3A, the DF linearization has two-order-of-magnitude higher accuracy than the first-order Taylor expansion, and can approximate the original full scale model with a reasonable root-mean-square error. This representation of nonlinear dynamics is markedly different from standard linearization methods. The DF linearization was also compared to the TPWL method. The figure shows an order-of-magnitude better result for the DF linearization compared to the sophisticated technique. Note that the TPWL does not yield a linear model since the state equation includes a product of two state functions. Therefore, superposition as applied with our DF linearization approach cannot be applied. As applied to the analysis of multi-cell ECM compaction, linear augmented equations describing single cell-ECM interactions were derived from DF linearization, and then converted to a reduced-order linear representation by transformation onto a basis of eigenvectors derived from simulated data set. Unlike model reduction of nonlinear dynamical systems, which still remains a challenging problem in the field [19–22], the model reduction of a linear system through DF Linearization is straightforward. It allows for the evolution of independent and auxiliary states to be described within a lower dimensional linear manifold. The resulting reduced order latent variable model is capable of reproducing nonlinear dynamics, and the linearized structure of individual models facilitated their integration to describe multi-cell behaviors. The prediction of collective behaviors of a group of cells was achieved by superposing contributions of individual cells represented by latent variables zc, k, which evolves based on their own dynamics in response to the global ECM state represented by latent variable ze. The linear representation of the collective multi-cell-ECM interactions manifests the two types of feedback actions by the individual cells. As shown in the block diagram in Fig 3A, the individual cells are exposed to the ECM forces represented by latent variable vector ze in two separate paths. The path through the cell polarity block and matrix B, leading to lamellipodia formation, can be viewed as an “active input” as addressed in [5]. This feedback path includes a cell’s internal decision as to which direction it extends lamellipodia. In contrast, the other feedback path through a gain matrix C does not have a high-level cell decision, but is reactive, playing a “passive role” [5]. These feedback interactions support the prior experimental work [5]. It is interesting to note that ECM compaction begins almost instantaneously, but the magnitude of compaction is rather limited. Once the “active” feedback loop is initiated in, the ECM compacts further, resulting in a large deformation. As the polarity dynamics are rather slow, the second stage ECM compaction does not start immediately. The time scale is determined by the constant κ involved in the polarity dynamics Eq 25. Using the proposed methodologies, we are able to reproduce intercellular mechanical interactions consistent with published experimental observations. In particular, the global compaction of gel volume via collective cell-contractile activities is characteristically different from local deformations of single isolated cells embedded within the same gel. Through study of emergent behaviors of groups of cells embedded in a 3-D ECM fiber network, we can advance our understanding of intercellular mechanical signaling during tissue formation [1–7]. There are a few limitations to our method, however. While the presented method can predict complex nonlinear behaviors, the method is still a type of approximation. Care must be taken with the validity period. In Fig 4C at the sample time of t = 50 minutes, the latent variable superposition simulation over predicts the volume shrinkage by 12%. With the current mathematical formulation, we have not yet incorporated the degradation of ECM fibers through matrix metalloproteinases. ECM degradation would be necessary to reproduce sustained movement and migration of the cells particularly in 3-D embedded matrices [32]. Since ECM degradation continuously changes the fiber connectivity through ECM remodeling, a methodology to update the node grid structure describing the ECM field would need to be developed. However, ECM degradation may not be necessary for predicting gel compaction since a cluster of cells remains stationary when contracting the surrounding gel [5]. Finally, in the current work, it was assumed that the cell’s polarity mechanism is a dominating internal response to mechanical cues. Cells change their internal state through a complex process of mechanotransduction and intracellular signaling. Incorporating these more complex mechanisms is an exciting avenue for future research. While the method has been developed and demonstrated for multi-cellular interactions with 3D ECM, the basic methodology is applicable to a broad range of systems where nonlinear dynamics of many interacting subsystems are prohibitively complex to compute.
Collective behaviors of multiple cells interacting through an ECM are prohibitively complex to predict with a mechanistic computational model due to its highly nonlinear dynamics and high dimensional space. We introduce a methodology where nonlinear dynamics of single cells are superposed to predict collective multi-cellular behaviors through a developed linearization method. We represent nonlinear single cell dynamics with linear state equations by augmenting the independent state variables with a set of auxiliary variables. We then transform the linear augmented state equations to a low-dimensional latent model and superpose the linear latent models of individual cells to predict collective behaviors that emerge from multi-cellular interactions. The method successfully reproduced experimental results of cell-induced ECM compaction.
cell physiology
focal adhesions
random variables
covariance
cell polarity
simulation and modeling
developmental biology
systems science
mathematics
molecular development
cellular structures and organelles
research and analysis methods
computer and information sciences
nonlinear systems
adhesion molecules
extracellular matrix
nonlinear dynamics
cell membranes
probability theory
cell biology
biology and life sciences
physical sciences
Multi-Cell ECM compaction is predictable via superposition of nonlinear cell dynamics linearized in augmented state space