A MULTIPHYSICS COMPUTATIONAL FRAMEWORK FOR UNDERSTANDING CELL AND MICROTISSUE MORPHOGENESIS

A sophisticated computational model is developed to consider different interactive parts between cells and its components and their local microenvironment. The present work is mainly focused on the modeling of the coupling of the stress and concentration of the signaling proteins within the cell domain. In this research, the fundamental aspects and details of a coupled Contraction-Reaction-Diffusion (CRD), is presented. The model accounts for diffusion of the proteins and mechanical equilibrium of the cells simultaneously while considering different subunits which are affecting the cell migration. For instance, cell-cell interaction, nucleus effects, focal adhesion distribution, anisotropic stress fiber formation, membrane tension, microtubule structure, and growth and retraction of the cells are considered. Collectively, because of the interaction of these different subunits, the cell works as a single migratory machine. The model fills the gap in coupled biomechanical and biochemical models for the biological cells and predicts both the instantaneous and the long-term dynamic behavior of the cells. In order to evaluate the proposed computational cell model, biological experiments such as cell migration, durotaxis, and collective cell migration has been simulated using the proposed computational model. The proposed model presents a simple mechanistic understanding of mechanosensing of substrate stiffness gradient at the cellular scale, which can be incorporated in more sophisticated mechanobiochemical models to address complex problems in mechanobiology and bioengineering. The proposed model and computer program is able to simulate long-term interaction of hundreds of cells with each other (e.g. cell-cell contact) and with the elastic substrate on a desktop workstation efficiently.

): a cell crawls from the stiffer side of the substrate toward the softer side and turned 90° at the interface (the dotted line is an approximation of the rigid-to-soft interface). (b) Traction stress under a circular cell crossing step-rigidity boundary (reprinted with permission from [14]) (c) Lamellipodia extension in a square cell (reprinted with permission from [15]). ..... 9     With computational simulations of cell migration and microtissue formation, optimal extracellular matrices can be designed to facilitate the recapitulation of micromorphologies of tissues at the cellular resolution in 3D bioprinting. Simulations are efficient tools for evaluating and optimizing the feasibility of specific designs.
It has been shown that mechanical and external forces play a key role in biological cells and their motility. Lo et al. reported the substrate rigidity-based cell migration, they discovered that fibroblast cells make a 90 degree turn at the boundary of the soft and rigid zones of the extra cellular matrix (ECM); this is referred to as durotaxis [1]. Barnhart et al. showed that the level of adhesion to substrate affects morphologies of the cells [2]. Kim et al. observed that a group of endothelial cells tend to fill the empty spaces of the substrate by exerting force and pulling out traction in the direction of empty spaces, which is called Kenotaxis [3]. All of the abovementioned phenomena are a few of many instances of the biomechanical behavior of 2 the biological cells. In general, mechanosensing is an interaction which is involved in a variety of biological processes such as cancer metastasis, wound healing, tissue formation, and embryological development.
Advancement of microscopy and cell imaging techniques including 2D and 3D force microscopy [4], multiple speckle microscopy [5] and monolayer stress microscopy [6] has enabled researchers to evaluate the mathematical models in cell mechanics considering the quantitative experimental data as the benchmark. In this context, developing a whole-cell mathematical model presents an ongoing challenge and has been the topic of intensive researches. For capturing the polarization of the cells during the cell migration (i.e., formation of the head and tail); inspiring by Turing's instability; some studies used the reaction diffusion of biochemical molecules [7], [8], [9], [10]. Although these models were able to mimic the polarization, they did not consider the invisible mechanosensing of the biological cells. Zaman et al.
considered the cells as a material point and were able to capture the motility behavior of a single cell in different stiffnesses of the ECM (i.e., higher stiffness associated with lower migration speed and lower stiffness associated with higher migration speed). However, their model cannot resolve the spatial distribution of the field variables such as traction stress, displacement and protrusion concentration within the single-cell domain [11]. Vernerey et al. proposed a static and rigidity sensitive cell model that can predict the stress fiber direction successfully [12]. A dynamic multi-cellular model is needed to consider the cell-cell interaction, effect of nucleus as well as the dynamic growth of the boundary at the front of the cell. Borau et al. presented a onedimensional model which focuses on the rigidity sensing of the cells and yet it is unable to consider the membrane tension, polarization, cell-cell interaction, and cell growth [13]. While each of the abovementioned models contributes to a better understanding of the cells behavior, lack of the studies on the biomechanics and biochemistry couplings of cell-microenvironment interactions has been identified, particularly computational models and methods that can describe and predict the dynamic process of microtissue formation.
Previously, our research was focused on developing a mechanical model for biohybrid cell contractility assays and studying the effect of thermal fluctuation on cell adhesion [14], [15]. In this study, an efficient and robust computational biomechanics model and software platform will be developed to simulate the dynamics of living cells at the whole-cell level. The proposed model establishes a firm link between protein distribution pattern and the traction stress in the cells. Specific biomechanics phenomenon including cell crawling and morphogenesis of cell monolayer tissues has been studied using the computational model.
The developed model integrates the biochemical and mechanical activities within individual cells spatiotemporally and it is mainly composed of four modules: mechanics of cytoskeleton, cell motility, cell-substrate interaction, and cell-cell interaction. In the cell membrane and cytosol domain reaction-diffusion equations of active and inactive diffusive molecules is formulated to model the protrusion and retraction protein concentrations. In the cytoskeleton domain elasticity equations is developed, and the mechanical stresses experienced by the cell has been solved. Finite element method (FEM) has been used to model the irregular shapes of cells and to 4 solve the resulting system of reaction-diffusion-elasticity equations. The weak coupling scheme between traction and protrusion and retraction protein concentrations has been adopted for this multiphysics problem. Automated mesh generation has been hired for re-meshing and to handle the element distortion in FEM due to the large shape changes of the cells.

Introduction
It has been shown that biological cells can sense and respond to a variety of mechanical cues of their microenvironment, such as matrix rigidity [1], matrix topology [2], matrix dimensionality [3], shear flow [4], interstitial flow [5], cell-cell and cell-matrix adhesions [6], and cell shape constraints [7]. These mechanical stimuli play a critical regulatory role in many biological functions such as cell proliferation 9 [8], cell motility [1], [7], and differentiation [9]. Understanding the mechanisms underlying mechanosensing has become the focus of intensive experimental and theoretical studies [10]- [13]. In the present study, we are interested in durotaxis, a termed coined by Lo et al. [1], which refers to the substrate rigidity-guided cell migration. They showed that (see Fig. 1a) when a fibroblast cell crawled from the stiffer side (i.e., the darker region) of the substrate toward the softer side (i.e., the brighter region), the cell made a 90-degree turn at the interface. The phenomenon of durotaxis (reprinted with permission from [1]): a cell crawls from the stiffer side of the substrate toward the softer side and turned 90° at the interface (the dotted line is an approximation of the rigid-to-soft interface). (b) Traction stress under a circular cell crossing step-rigidity boundary (reprinted with permission from [14]) (c) Lamellipodia extension in a square cell (reprinted with permission from [15]).
Because the importance of durotaxis in physiology and pathology, the molecular and subcellular mechanisms underlying mechanotransduction has attracted considerable attention [16]- [19]. Biomechanics models where single cells and cellsubstrate linkages were modeled as elastic springs or elasticity theory, have been developed to account for rigidity sensing [13], [20], [21]. These models were only applied to the scenario where the individual cell was treated either as a point mass [20] or an small volume element [21]. In a series of studies [22]- [24], by modeling the individual cells or stress fibers as force dipoles distributed in continuum elastic 2D or 3D substrates, the researchers developed biomechanics models to interpretate mechanosensing mechanisms and to study the effect of mechanosensing on cell shape, stress fiber orientation, and synchronized beating of cardiomyocytes. On the other hand, the effect of durotaxis on cell migration dynamics on the long time have also been studied. For example, in the cell migration model by [25], substrate rigiditydependence is taken into account by assuming focal adhesions are correlated with substrate stiffness. In the single cell migration model by [26], substrate rigiditydependence is considered by assuming differential cell-substrate adhesions strengh on substrates with different rigidity. Thus, in these durotaxis models, substrate-rigiditydependence is used as the assumption in the cell migration models. To the best of our knowledge, in previous models, the spatial distribution of the substrate rigidity within the single cell domain is constant. The substrate rigidity is either changed for the whole cell or changed only when the cell moves from one location to another (during migration). The main difference between these models and our model is that we examine how a single cell senses the local substrate rigidity difference within the single cell domain. Therefore, our model can provide a more direct interpretation on how the cell sense the rigidity gradient.
In their original paper [1], Lo et al. proposed a two-part theory for the detection of the spatial gradient of the substrate rigidity as follows. In the first part of their theory, the cytoskeleton-focal adhesion-substrate linkages are considered as elastic springs (with spring constant k for the same amount of elastic energy input U from the active actomyosin contraction) to pull these springs, the spring force F enerated is larger at the stiffer region of the substrate underneath the cell (Because thus for the same U larger k results in larger F In the second part the theory, the stronger force leads to a higher level of activation of force-sensitive proteins through conformational changes, which in turn leads to migration-related cellular responses such as upregulation of lamellipodia extension. The second part of the theory is referred to as mechanotransduction [27] in the literature. This two-part theory is directly supported by other experimental observations. For example, in a work by Breckenridge et al. [14], traction stress under a circular cell crossing step-rigidity boundary were measured using elastomeric micropost arrays. They found that the traction stress is higher on the stiffer half of the circular island (see Fig. 1b), which supports the first part of the theory. In another work [15], the authors found lamellipodia grow preferentially from the corners of square cells (see Fig. 1c). The corners of convex polygonal shapes are known to be the spots where high traction stress is generated when the cell contracts [28]. Together these findings support the second part of the durotaxis thoery that larger focal adhesion force leads to more lamellipodia extensions.
In this work, we are concerned with the first part of the theory: how exactly is the larger focal adhesion force generated in the part of the cell adhering to the stiffer region of the substrate? The assumption of Lo et al. that the same energy is provided for pulling is not without pitfalls, since it is not straightforward as how the generation of the same mechanical energy is ensured at the different sub-regions of the cell by the cell's active contractile apparatus. Another (and easier) approach to calculate the force is to consider the static equilibrium of the cell. The migrating speed of fibroblast cells is very slow (~1 µm/min), considering the stress fibers are in a state of isometric tension, thus the cell at any time instant can be considered to be in a quasi-static equilibrium. Therefore, the static equilibrium holds for the whole mechanical system composed of the cell and the elastic substrate.
Using the method of static equilibrium, a simple generic model based on active matter theory has been devised by Marcq et al. [29], in which the cytoskeleton was modeled as two parallel elements (one passive spring and one active contractile element), and the 1D cell connects to the substrate springs only at the two ends (see Fig. 2a). Their model is sufficient to explain the experimental findings [17], [30] where the magnitude of the traction stress increases with the substrate rigidity.
However, it cannot explain the rigidity gradient sensing (i.e., different traction stress at the two ends of the same cell): the static equilibrium implies that the adhesion forces 13 at the two ends of the cell should be the same, regardless of disparate substrate-spring stiffness. This is the paradox that was raised in a review by [31].
The assumption that the 1D cell only adheres to the substrate at the two ends oversimplifies the problem. In fact, by dropping this assumption, the abovementioned paradox can be resolved. Considering that the cell adheres to the substrate in the whole cell domain, in the present study, we show that the static equilibrium of the cell is sufficient to yield the rigidity gradient-dependent traction force distribution. The remainder of the paper is organized as follows. We first present the simple elasticity model for 1D and 2D cell adhering to an elastic substrate. For the 1D cell, we will derive analytical solutions and present results from the parametric studies. For the 2D cell, we will use the finite element method (FEM) to numerically solve the equilibrium equations. We then compare the modeling results with the three experimental observations listed in Fig. 1.

1D Model
We first present a 1D model of a cell adhering to an elastic substrate. As shown in Because the substrate is modeled as isolated springs (i.e., elastic interaction within the substrate is neglected), the substrate considered in our model can be thought as the elastomeric micropost arrays [14], rather than a conitnouum elastic substrate. The active actomyosin contraction shortens the 1D cell and the shortening is resisted by the passive compoment of the cytoskeleton and the substrate (see the schematic in Fig. 2b).
Constitutive relations of the cytoskeleton have been previously studied intensively. Time-dependent constitutive relations based on Hill's law of muscle contraction have been devised previously to capture the dynamic process of actomyosin contraction, such as for stress fibers [13], [17], [32] or for myofibrils [33], [34]. In this work, for simplicity, the final state of the dynamic models when contraction stress reaches isometric tension and strain rate becomes zero is considered, which yields a time-independent consituttive relation for the 1D cytoskeleton: where is the overall cytoskeleton stress, is the Young's modulus of the passive component of the cytoskeleton, is the isometric tension due to the active actomyosin contraction, is the strain, is the displacement along the axis of the 1D cell (i.e., -axis).

16
The static equilibirum equation of the 1D cell is (2) where is the thickness of the cell and is assumed to be a constant for simplicity, is the traction stress exerted on the substrate by the cell. Because the focal adhesion is connected to the substrate spring in series, is also the stress experienced by the focal adhesion. Therefore, we here use the phrases "focal adhesion stress" and "traction stress" interchangeably in this paper. Traction stress x T an be calculated as (3) where is the equivalent spring constant of the cell-substrate linkage composed of the focal adhesion spring and the substrate spring, as illustrated in Fig. 2b. Because the two springs are in series, we have To model the rigidity gradient, we define step changes in substrate rigidity by (5) where  efines the ratio of rigidities of the two regions, which can be regarded as the gradient strength. Without loss of generality, the left half (i.e, 0 x  is considered to be softer than the right half (i.e., 0 x  Therefore, 1   is imposed in our parametric studies. The stress-free boundary condition applies at the two ends:

2D Model
To apply the model to cells cultured on 2D surface of elastic substrate, we extend the 1D model to the 2D. For the 2D model, Eq. (1)-(3) become where the indicial notation is used, summation over repeated indices is adopted, , , and are stress tensor, deviatoric strain tensor, strain tensor, respectively, and are shear and bulk moduli of the cytoskeleton. the substrate rigidity gradient is modeled by defining as a function of Cartesian coordinates . The finite element method (FEM) [35] is used to numerically solve the differential equations of the 2D model, where 3-node triangle element is used for the spatial discretization.

1D cell, when α = 1
When the gradient strength parameter 1   , meaning uniform rigidity underneath the cell, the 1D model can be readily solved for the displacement x u and traction x T as follows. Using the strain-displacement and constitutive relations, Eq. (2) can be rewritten as an linear second-order differential equation: , 0 where cs k is given in Eq. (4). Figure 3a plots END T as a function of ECM k which shows that the traction stress reaches a plateau when ECM k   which implies there is a saturation value of traction stress or force at large substrate rigidity. This result has been previously shown in experiments and models [17], [29], [30] Mathematically, this is simply because Mechanically, this is because two springs in series is softer than any of the two springs. Therefore, when the substrate becomes rigid, the spring stiffness of the cell-substrate linkage becomes equal to the focal adhesion spring.

1D cell, when α >1
When  >1, a step change of rigidity is present underneath the single cell (the left half is always softer than the right half). The analytical solution can be derived similar to the case of  =1. Figure 3b and 3c show the solutions of traction stress and cytoskeletal stress , respectively. For the traction stress, positive sign means rightward pulling and negative sign means leftward pulling. Clearly, the traction stress magnitude is maximal at the cell edge on the stiff region (i.e., at the position ). With continuous adhesion to the substrate in the whole cell domain, traction stress is redistributed so that on the stiffer side the traction stress is within a shorter range but higher magnitude on average. On average, the higher cytoskeleton stress is generated in the stiff region compared to the soft region. Note that traction stress is discontinuous at the interface (i.e., 0 x  ). This is simply because the substrate rigidity is assumed to be discontinuous at the interface in our model (see Eq. 5, there is a step change of rigidity across 0 x  ). If we assumed a linear-varying rigidity gradient, the traction stress would be continuous. Parametric studies were conducted to ascertain the sensitivity of the modeling results to the parameter values. We define the difference between the traction stresses at the left and right ends of the 1D cell as denotes the absolute value. Quantity Δ x T represents the difference between the traction stresses on the soft and stiff regions of the substrate. Figure 4 plots the Δ x T as a function of  for different values of E and s k In both Fig. 4a and 4b, we see that the traction stress difference increases with  which implies that the gradient strength play an important role in durotaxis [36], [37]. In Fig. 4a, we can see that Δ x T increases as E decreases, meaning softer the passive component of the cytoskeleton, larger difference of traction stress is produced. When the passive component of the 21 cytoskeleton becomes stiffer, less contractile force is transmitted to the focal adhesion and consequently weaker dependence of focal adhesion stress on the substrate rigidity.
In Fig. 4b, we can see that the ratio between the relative difference of traction stress between the soft and stiff regions (i.e., , where = )increases with decreasing substrate rigidity s k which implies that if the mechanotransduction process detects the relative difference of traction stress, then softer substrate promotes durotaxis.

Computational results for 2D cells
First, we show in Fig. 5a the FEM simulation results from the 2D cell model for a circular cell crossing a step-rigidity boundary (i.e., the upper half of the cell adheres to soft micropost arrays, the lower half of the cell adheres to stiff micropost arrays). As shown in Fig. 5a, the traction stress is higher on the perimeter of the cell, and it is 22 higher on the lower half (stiff substrate) compared to the upper half (soft substrate).
The displacement is slightly higher on the upper half than the lower half. Our modeling results in Fig. 5a are in good agreement with in the experimental results [14] shown in Fig. 1b (if we neglect the random noise in the experiment). Therefore, both the experiment and our model show that the cytoskeleton contraction in the single cell generates higher traction stress on the stiffer region of the substrate underneath the cell.
Second, we show in Fig. 5b the model prediction of traction stress for the square cell. The traction stress concentrates to the edge of the square cell, and maximizes at the corners. This modeling result is correlated with the experimental data by [7], [15] shown in Fig. 1c where the lamellipodia extensions were localized to the corners of square-shaped cells. This correlation supports the durotaxis theory proposed by [1]: larger focal adhesion forces at the corners of the square cell are converted into protrusion signals via molecular mechanisms of mechanotransduction, which eventually lead to stronger lamellipodia extension. In the case of durotaxis, higher traction stresses are in the rigid side of the substrate and essentially the protrusion signals will be amplified in the rigid side rather than in the soft side.  Fig. 1a to be the interface between the soft and stiff sides of the substrate. Figure 5c shows the traction stress distributions corresponding to the experimental images of Fig. 1a at the different time instants. One can see that the traction stress for the lamellipodia on the rigid side (solid arrowhead) is larger than that of the lamellipodia on the soft side (hollow arrowhead). If the larger force is converted into more protrusion signal, the lamellipodia on the right side (solid arrowhead) will become the dominant one, which eventually leads to the turning of the cell at the step-rigidity boundary. Note that the highest traction stress spot at the tail of the cell at the beginning (see Fig. 5c) does not result in a leading head is probably because the memory (in the molecular constitutes) of the head-to-tail polarization [38], i.e., the new head will most likely to form near the original head.

Conclusions
In this work, we use a simple elasticity mechanics model to predict the traction stress (i.e., focal adhesion stress) for single adherent cells on the elastic substrate with rigidity gradient. The model predicts larger traction stress (i.e., larger focal adhesion stress because the traction stress is equal to the force experienced by the focal adhesion) on stiffer region of the substrate underneath a single cell. This minimal mechanics model provides a plausible answer for the first part of the durotaxis theory proposed by Lo et al. [1]: how exactly is the larger focal adhesion force generated in the part of the cell adhering to the stiffer region of the substrate? We found that the principle of static equilibrium alone provides a mechanistic explanation to this question. We think our model has resolved the paradox that was raised in a review by [31], which states that a static model cannot explain the rigidity sensing of a cell. Our model can be incorporated in more sophisticated mechanobiochemical models to address complex problems in mechanobiology and bioengineering.

Introduction
Cell migration, an intriguing phenomenon involved in tissue formation and cancer metastasis, has long been the subject of intensive investigation in the fields of 35 biophysics and cell biology [1]. The machineries that drive migration, and the signaling networks that control the migration machineries have been studied intensively [2]. In vitro 2D cell migration experiments revealed that cells of different types, such as fibroblasts, keratocytes, and neurons, exhibit different migration/spreading behaviors [3]. On the other hand, different types of cells share some fundamental characteristics. In general, single cell migration can be described as a coordinated and integrated process of different modules ( Fig. 1): cell polarization (i.e., front-and-rear formation), protrusion of lamellipodia/filopodia/lobopodia, invadopodia formation and proteolysis of surrounding ECM, formation of new adhesions in the front, releasing of aging adhesions at the rear, and cytoskeleton/membrane skeleton contraction to move the rear forward [2], [4].
Inspired by Turing's reaction-diffusion model of diffusive biochemical molecules, mathematical models were developed to study cell polarization [5], [6], and cell morphogenesis in migration [3], [7]. Particularly, the reaction and diffusion of intracellular signaling molecules have been interpreted as to form networks [8], [9].
Cells with this internal network system are able to spontaneously polarize and make persistent random walks in the absence of external cues [3], and to carry out directed movement when biased by external signal gradients (i.e., chemotaxis).
On the other hand, biomechanics has been shown to play a critical role in many biological functions such as cell motility [10]- [12]. It has been postulated that the mechanosensitive proteins change their conformations when stretched by mechanical forces, and the conformational changes open up hidden active sites for binding with other molecules, which in turn results in specific chemical reactions. For example, 36 stretching of talin protein resulting in the increased binding of vinculin to talin [13], which may be a molecular mechanism underlying the force-dependent focal adhesion maturation [14], [15], [16]. Similarly, stretching of α-catenin has been shown to induce enhanced vinculin binding in cell-cell adhesion [17], [18]. Mechanical forces also regulate the actomyosin stress fiber formation [19]. Therefore, mechanical and geometrical properties of cells and their microenvironments are not merely passive consequence of biochemistry. In steady, they are important regulators of biological processes such as cell migration.
Biochemical models based on reaction-diffusion equations lack the consideration of mechanotransduction thus cannot capture mechanosensing in cell migration.
Biomechanics models lack consideration of biochemical signaling and thus fail to account for biochemical processes. A thorough understanding of cell migration will require the elucidation of how the mechanical and biochemical events are spatiotemporally integrated at the cellular scale. In this work, we develop a contraction-reaction-diffusion model for cell migration by integrating the mechanical force generation and reaction-diffusion of biochemical molecules at the whole-cell scale. Our overarching hypothesis is the following: the biomechanical and biochemical signals form mechanobiochemical feedback loops. Through the mechanobiochemical feedback loop, cell migration and cell shape can be regulated by a variety of mechanical cues, such as cell shape [20] and substrate rigidity gradient [21]. asymmetry, protrusion and retraction of cell body, and cell-substrate interaction, which will be described below.

Cytoskeleton module
The elasticity model of the cytoskeleton A rather simple mechanics model of cytoskeleton is adopted here. Because the migration of biological cells in their solid or fluid microenvironment is at the low Reynolds number [22], the inertia force can be neglected. At each time instant, the cell can be considered in a quasi-static equilibrium. The equilibrium equation of the 39 cytoskeleton in the domain are written by using Cartesian tensor notation (summation over repeated indices is adopted hereinafter) as (1) where is the Cauchy stress tensor of the cytoskeleton ( and takes values of 1 and 2 in 2D), is the thickness of the cell, which is assumed to be uniform throughout the cell. Denoting the area of the cell, the volume of the cell is calculated as (2) The cell volume is assumed to be conserved in the present study, so where is the strain tensor, is the displacement, is the deviatoric strain tensor, and are shear and bulk modulus of the passive network, is the active isometric tensile stress (ITS) tensor from the active part of the cytoskeleton, which will be defined later in Eq. (6). Use of the smallstrain Hooke's law in Eq. (3) for the large deformation that occurs during cell motility deserves some explanations here. In this model, when solving the elasticity problem for a migrating cell, at each time instant we treat the current configuration as the 40 stress-free state, and displacement and strain are measured from the current stress-stress configuration, rather than the reference state at an initial time of the whole migration process. Note that the displacement and strain concepts here are different than those in the conventional finite-strain elasticity theory. This is an ad hoc treatment based on the assumption that the dynamic bonds that forms the passive network of the cytoskeleton remodels fast enough to release the passive stress in the cytoskeleton. The traction stress is assumed to be linearly proportional to the "instantaneous" displacement of the cell (4) where is the spring constant of the cell-substrate linkage that will be defined later in Eq. (18).

Stress-fiber structure tensor
To account for the anisotropic fiber formation, a previously defined mathematical model for myofibril orientation in cardiomyocytes is used [20]. A second-order tensor , referred to as the stress-fiber structure tensor, is defined and its time evolution is described by (5) where and are the stress fiber activation and deactivation rates, respectively. Here, is a model parameter, is the cytoskeletal tension, denotes Heaviside function and is defined as: =1 when and =0 when . Denoting the maximal eigenvalue and the corresponding eigenvector of by and , respectively, the ITS tensor is defined as 41 (6) where , , and denotes the baseline, retraction signal-associated, and stress fiber-associated contractility, is the concentration of retraction signal that will be introduced below. The dyadic product of unit vector produces the tensor , which has its only non-zero-principle-value principle direction along . This anisotropic cytoskeleton model has been used to explain the pattern formation of myofibrils in single cardiomyocytes [19].

Cell motility module
In this work, we adopt a similar modeling concept as Satulovsky et al. [3] where a few phenomenological variables are used to represent the concentration of various proteins involved in cell migration.

Reaction-diffusion of protrusion and retraction signals
Previous studies indicated that the active forms of protrusion and retraction signals are membrane-bound proteins [7], [23]. Two phenomenological variables and are defined in the physical domain of the membrane to account for the area concentration of active form of protrusion (e.g., Rac, Cdc42) and retraction (e.g., ROCK) signals, respectively [3]. Here variables and are normalized by their saturation values respectively thus are in units of µm -2 and with the maximal value of 1 µm -2 . Their time evolution equations are defined as To be simple, and are assumed to be uniform in the cytosol and are calculated by the following two equations, respectively, where and are model parameters denoting the total volume concentrations of both active and inactive forms for protrusion and retraction signals, respectively, and are the spatial mean values of and , respectively. In Eq. (9), multiplying by cell thickness to convert a volume concentration to an area concentration is based 43 on the assumption that the cell is flat and the diffusion in the cell thickness direction is instantaneous.

Movement of the cell (i.e., protrusion and retraction)
The movement of the cell consists of the cell protrusion caused by the actin polymerization at the leading edge of the lamellipodia and the passive retraction as a result of active cytoskeleton contraction. A protrusion velocity is defined as a function of the protrusion signal at the cell edge as (11) where model parameter represent the intrinsic level of membrane protrusion speed, parameter is the threshold of protrusion signal above which membrane protrusion occurs, parameter sets an upper limit of the cell area, is the outward normal unit vector at the cell edge. The retraction velocity is assumed to be proportional to the instantaneous displacement of the cytoskeleton as, where c r is a model parameter.

Cytoskeletal asymmetry
Experimental studies have revealed that the cell polarity (i.e., head-and-tail pattern) is maintained through the long-lived cytoskeletal asymmetries including microtubules [24]. To incorporate the cytoskeletal asymmetries in the model, a vector is defined to represent the polarity of the asymmetric cytoskeleton and its time evolution equation is defined as, 44 (13) where is a model parameter, vector is a vector defined based on the protrusion signal, Where is a model parameter, is a unit vector and is the position vector of the edge points relative to the center of the cell, is the length of , is the differential angle corresponding to the differential arc length, where is the angle coordinate of the edge point in the polar coordinate system with the cell center as the origin. As implied by Eq. (13), in the steady-state ( =0), the cytoskeleton asymmetry vector is equal to the vector . The cytoskeletonasymmetry function in Eq. (7) is defined with the angle of the vector as, where is a model parameter.

Cell-substrate interaction module
To incorporate the dynamic remodeling of focal adhesion, a phenomenological variable is defined to describe the density distribution of focal adhesion-associated proteins (e.g., integrins, talins, vinculins, etc.). Variable is normalized by the saturation value, thus ranges from zero (no integrin-mediated cell-substrate adhesion) to one (mature FAs). The time evolution of is described by , where , , , and are the rate constants for the spontaneous, autoactivation, protrusion signal-dependent, and stress-mediated focal adhesion formation, respectively, is a decay constant, denotes the magnitude of the traction stress, and are model parameters, and is the strength of random noise. Here the redistribution (e.g., via active transportation and passive diffusion) of unbound focal adhesion proteins is assumed to be faster than other time scale of focal adhesion formation, the unbound focal adhesion protein density in the membrane is simply computed as (17) where is the mean value of in the membrane domain, and represents the average density of the total amount of bound and unbound focal adhesion proteins.
Denoting and the equivalent spring constants of the focal adhesion and the substrate, respectively, the spring constant of the cytoskeleton-substrate linkage is given as (18) The mechanics of the cell is coupled to the dynamics of focal adhesion remodeling through the spring stiffness by the following relation (19) where is the maximal stiffness when the focal adhesion density = 1 µm -2 (i.e., mature focal adhesion).

Numerical implementation of the model
The cell monolayer model is implemented in an in-house code using the finite element method, where Lagrangian mesh is adopted and 3-node triangle element is 46 used. In the simulations of the movement of the cell, the nodal spatial coordinates are updated based on Eq. (11) and Eq. (12). An auto mesh-generating algorithm based on Delaunay triangulation is utilized to perform re-meshing when mesh distortion occurs.
Mesh transfer for the field variables is performed between the old and new mesh.

Mechanosensing at the whole-cell scale
As illustrated in Fig. 2

Establishment of polarity of the cell with the reaction-diffusion submodel
Cell polarization (i.e., forming head and tail) is critical in cell migration to achieve directed movement. The spontaneous polarization has been thought as a pattern formation in reaction-diffusion systems [3], [23], [25]. We here define the system of equations consisting of Eq. (6)- (9), where and are set to be zero, as the reaction-diffusion submodel. In this reaction-diffusion submodel, which were previously proposed by Maree et al [7], the protrusion signal and retraction signal inhibit each other through the and terms, respectively. As studied by Maree et 48 al, this simple mutual-inhibition model can induce spontaneous polarization. Figure 3 shows the simulation result of the reaction-diffusion submodel. As shown in Fig. 3A (top row), starting with a randomly perturbed initial state, the protrusion signal in a circular cell spontaneously polarizes, i.e., spatially separates into two zones: high and low regions. Because of the auto-inhibition, retraction signal distribution also polarizes.
49 As previously showed by Maree et al. [23], the cell shape (i.e., the shape of the mathematical domain of the reaction-diffusion equations) has an important effect on the spatial patterns formed. They concluded that at the steady state, the length of the interface that separates the high and low regions is minimized. Our simulation results agree with their conclusion. As shown in Figure 3B, the interface in the elliptic cell is initially setup to be parallel to the longer axis of the ellipse (i.e., the initial distribution of the protrusion signal is a gradient from high in the left to low in the right). Over time, the interface rotates and eventually aligns with the shorter axis of the ellipse. To see how various model parameters, influence the steady-state protrusion and retraction distributions in the reaction-diffusion submodel, a rectangular (with an aspect ratio of 1:3) cell shape is used to simulate approximately a quasi-1D domain. The simulation results for parameter , , , and are plotted in Fig. 3C-F. One can see that change of these parameters can all shift the position of the interface and the peak value of at the steady state.

Focal adhesion stress-dependent protrusion signal distribution
The reaction-diffusion submodel described above can account for the polarization of the protrusion signal, but it cannot explain the phenomenon observed in the previous study [15] in which the membrane protrusion localized to the four corners of the square cell. As shown in Fig. 4A, the reaction-diffusion submodel alone predicts a polarized pattern for the protrusion signal distribution in a square cell.
Our previous studies showed that localization of the traction stress (which is equal to the focal adhesion stress) at the corners of the square cell is simply due to the mechanics principle of static equilibrium of an elastic body [19]. Based on that, we argue that protrusion signal and focal adhesion assembly can be enhanced by the mechanical stress in the focal adhesion. This hypothesis is implemented in our model by introducing the term in Eq. (6) and the term in Eq. (16). (16)), larger leads to bigger (Eq. (18) and (19)), bigger results in larger traction stress (Eq. (4)).

The role of the cytoskeleton asymmetry to the persistent migration
In this model, we introduce a cytoskeleton-asymmetry function to be able to explicitly control the directional persistence of migration. Figure   adhesion were manually tuned such that both the pattern and magnitude of the predicted traction stress distribution best match the experimental results (see Fig. 6).
As a result, we found that one set of model parameters can be found to yield good matching for the most of the time frames in the whole course of migration period presented in Fig. 6.

Simulation of durotaxis
Durotaxis is a term coined by Lo et al. [10], which refers to the substrate rigidityguided cell migration. They showed in the in vitro experiment where a fibroblast cell crawls from the stiffer side of the substrate toward the softer side, the cell made a 90degree turn at the interface to avoid migrating into the softer region. A conceptual two-step theory consisting of the force generation and mechanotransduction has been proposed previously by Lo et al. to explain the durotaxis. A simple mechanics model has been presented by us previously to explain how exactly the larger focal adhesion stress is generated at the stiffer region of the substrate [14]. We showed that static equilibrium of the adherent cell along can yield the disparate traction stress on regions of different rigidity. In this study, we integrate the elasticity model with the reactiondiffusion equations to form a contraction-reaction-diffusion system.
Here the dynamic model of single-cell migration is used to reproduce durotaxis phenomenon in silico. Cells started as a polarized circular shape and placed on the stiffer region of the substrate. The cell then crawls toward the softer region (i.e., left side). The results from two simulations are presented in Fig. 7. The only difference between these two simulations is the stress-dependent protrusion signal parameter : relatively high for the simulation I (top three rows) and low for the simulation II (bottom three rows). Note that parameter regulates the level of force-dependent activation of protrusion signal. At relatively low value of (top three rows), the cell crosses the interface with a slight change of cell shape. At high value of , the cell makes a turn at the stiff-to-soft interface and then crawls along the interface. The turning of the cell at the interface is caused by the positive feedback loop mentioned previously in Fig. 4: larger substrate stiffness leads to larger traction stress, larger traction stress leads to higher level of focal adhesion and protrusion signal, which results in change of migrating direction at the interface. These simulation results show that our model can successfully simulate durotaxis phenomenon.
Our simulations demonstrate that the cell can sense the non-specific mechanical cues of its microenvironment through a mechanobiochemical system (see Fig. 2). The key and the starting part of this system is the active contraction of the cell. Without the actomyosin contraction, no forces will be generated and the mechanosensors will not be activated.

Cell shape
Cell shape is an emergent property of a cell during cell spreading and migration.
Here we demonstrate how the cell shape can be changed by varying the model parameters. The parameter space of the model was searched to identify the cells of different shapes. To quantitatively characterize the cell shape, we define two dimensionless numbers: the roundness and branchness numbers. The roundness number, denoted by , was defined as where is the area of the cell, is the radius of the circle that circumscribe the cell boundary [3]. The roundness number takes maximal value of 1 when the cell shape is a perfect circle and is smaller than 1 for any other shapes. It is useful in distinguish between elongated and rounded shapes, but may fail to distinguish between the elongated and dendritic shapes. The branchness number is defined as , where is the perimeter of the cell. The branchness number has a lower bound of when the cell shape approaches a strip with zero width, and becomes large when the cell shape is dendritic.
Both the brute-force search and genetic algorithm were used for the parametric study. As shown in Table 1, three characteristic shapes were identified in the parameter space search: elongated, rounded, and dendritic shapes, with their roundness and branchness numbers listed. The values of some key model parameters that corresponds to these three characteristic shapes are also listed in the table. Note that the parameter values are given as ranges, indicating the regions of the parameter space where these characteristic shapes emerge.  One of the most common observation from analyzing the simulation results was that the total focal adhesion is an important parameter and plays crucial role in cell morphology. Altering the focal adhesion, changes the cell shape drastically.
In this work, we developed a computational model that integrates the biomechanics and biochemistry of the cell spatiotemporally at the whole-cell scale.
Biomechanical and biochemical events and processes were treated as modules, between which cross-talks are defined. We have shown that the reaction-diffusion submodel can simulate cell polarization (head-to-tail formation), the contraction-  Email: hongyan_yuan@uri.edu

Introduction
Tissue/organ morphogenesis is a complex process occurring at multiple scales.
Focusing on the whole organ scale, considerable research has been devoted to elucidation of the physical principles underlying the formation of the overall morphologies of organs [1]- [3], as well as the nutrient consumption and transport in 69 bioreactors [4] for tissue engineering. In these whole-organ level studies, information at the individual cell level has been homogenized or ignored. At the other extreme of the length scale, the genetic and molecular causes that dictate the tissue/organ formation have been intensively studied [5]. There is gap between our understanding of how phenotypic morphologies at the organ level emerge from genetic information.
Studies at the cellular and microtissue level play an indispensable role to bridge these two scales.
The phenotypic morphologies of cells including cell shape and cytoskeleton architecture, cell-ECM and cell-cell adhesions, can be best seen by comparing four types of tissues: muscle tissue, nerve tissue, epithelial tissue, and connective tissue.
Each of these different tissues exhibit characteristic morphologies in cell shape and cytoskeleton architecture. These four basic types of tissues are arranged spatially in various patterns (e.g., sheets, tubes, layers, bundles) to form organs. Gene expression only dictates what proteins to make and subsequently what biochemical reactions to carry out, the emergence of spatial morphologies must be determined by biomechanical principles and the coupling between biomechanics and biochemistry [6]- [8]. Mechanobiochemical coupling is exemplified by the recent discoveries in the field of mechanobiology. Cellular functions including cell migration and cytoskeletal dynamics that are closely related to cell morphogenesis, have been shown to be regulated by various mechanical cues such as matrix elasticity [9], matrix topology [10]- [15], matrix dimensionality [16]- [20], cell-ECM/cell-cell adhesions [21], and cell shape constraints [22]- [25]. Therefore, mechanical and geometric properties of 70 cells and their microenvironments at the length scale comparable to single cells can have a dominant effect on the microscopic tissue morphology.
Mathematical models based on reaction-diffusion equations at the cellular scale were developed to understand spatial pattern formation in the context of cell migration, such as cell polarization [26], [27] and cell morphogenesis [28], [29].
However, biochemical models lack the consideration of mechanotransduction thus cannot adequately capture cell morphogenesis. Biomechanics models were developed to interpret specific aspects of cell spreading, for example, the distribution patterns of traction force [30], cell adhesion [31], and cytoskeleton dynamics [32]- [36]. In contrast, biomechanics models lack consideration of biochemical signaling and thus fail to account for biochemical regulations. A thorough understanding of cell and microtissue morphogenesis will require the elucidation of how the mechanical and biochemical events are spatiotemporally integrated at the cellular scale.
In this work. We extended the single-cell model (See Chapter 3) to multicellular monolayer model by adding a module of the cell-cell interaction. Finite element method is used to solve the resulting system of partial differential equations and the model was implemented in an in-house MATLAB code.

Nucleus deformation and movement
The nucleus is modelled as an elastic structure that deforms upon the compression of the cell membrane and moves with the cytoskeleton. In the present model, there are no mechanotransduction associated with the nucleus. Rather, the nucleus is a passive material and can resist deformation and contribute to the shape of the cell in cases cell are elongated or compressed. In the finite element-based numerical implementation, the nucleus is discretized into networks of nonlinear springs connected at the nodes.
The configuration of the network is updated using Newton's equation of motion of the nodes. The numerical implementation concerning the passive nucleus model will be 73 presented in a later publication where the numerical algorithm and software package is described in detail.

Cell-cell interaction module
In cell monolayers where cells are connected mechanically by cell-cell adhesion, the static equilibrium of the cells depends on the cell-cell contact [37]- [40] in addition to cell-substrate adhesions. To simulate the dynamic process of formation and dissociation of cell-cell adhesion, a stochastic model is used to determine the binary state of the cell-cell adhesion as follows. When cell-A and cell-B is in close contact, the state of the cell-cell adhesion can be either "on" or "off". The "on" state indicates that the cell-cell adhesion is established. The "off" state indicates that although two cells are in close contact, they do not adhere to each other. The probability of the "off" state per unit edge length and unit time is denoted by cell-cell break rate parameter .
when the cell-cell adhesion state is "on", the stress between cell-A and cell-B, , is calculated as (1) where and are the displacement of cell-A and cell-B at their edges (where the cell-cell adhesion is formed), respectively, and is the spring constant of the cell-cell linkage.

Mechanobiochemical coupling
These different modules are coupled through the mechanics of the tissue. As illustrated in Figure. 1B, through molecular scale mechanotransduction pathways, 74 mechanical stresses in the cell are converted into biochemical activities, which in turn regulate the assembly/disassembly of macromolecular of the cell. The macromolecular assembly and disassembly alter the structural, geometrical, and material properties of the cell, which, according to the continuum/structural mechanics theory, will subsequently change both the internal stress (cytoskeleton stress) and stress at the boundary (cell-matrix and cell-cell adhesion stress). Thus, mechanics of the cell, biochemical activities, and macromolecular assemblies are coupled through mechanobiochemical feedback loops as depicted by the arrows in Figure. 1B.

Numerical Implementation
The cell monolayer model has been implemented in an in-house code. Finite element method has been used to solve partial differential equations resulted from reaction-diffusion and elasticity equations. Explicit Euler scheme has been used for the time integration. Due to extremely large deformation experienced by the cells during the cell migration, Lagrangian mesh has been adopted and 3-node triangle element has been used.
In terms of the algorithm implementation, we have 4 major part; data input, main function, cell migration, remeshing. Figure 3 shows a detailed algorithm that being used for the proposed model. The flowchart shows the main flow of the algorithm on the left side and on the right side it shows the main tasks within every loop. The flowchart introduces the core operations and almost each core operation itself consists of several subroutines (there is a list of major subroutines in the Appendix A). We 75 also use a flags-based coding, meaning different modules of the code can be either turned on /off or switched based on the scenarios. An example for the former is turning the cell migration on or off and an example for the latter is switching between initial distribution of the signals from random to constant or gradient pattern.
First, a data structure name para is used for the parameters input (see table 1 in appendix B for the parameters used in this study), the data structure will be initialized for feeding into the main function. Note that for simulating different scenarios (e.g.,  78

Case study: Collective cell migration in monolayer tissue
Collective cell migration has been studied in vitro where a confluent monolayer of cells crawls on a flat 2d substrate [41], [42]. Here we conduct the in-silico modeling of cells crawling in monolayers, as shown in Fig. 3. Total of 26 cells are confined in an adhesive region of a circular shape with a hole in the middle. The inner and outer radius of the adhesive region are 30 um and 83 um, respectively. To study the role of intercellular adhesion in collective cell migration, we performed two simulations: case-I (cell-cell adhesion is turned off) and case-II (cell-cell adhesion is on). The dynamic simulations start with circular cells seeded onto the adhesion region.

Conclusions
In this work, we developed a 2D computational model for studying cell morphogenesis in monolayer tissues. Because of the complex nature of the living cell, fidelity. This computational model will be used to elucidate and understand the morphological pattern formations in the microtissues that consist of many cells.

Parameter space search
The parameter space and mathematical formulations will be searched to identify the sub-spaces in which the computational model predicts characteristic behaviors of each type of differentiated cells, as well as the characteristic microscale morphologies of each type of tissues (e.g., bundling in muscle tissue, branching in nerve tissue, polarization in epithelial tissue).

90
The genetic algorithm will be used to estimate the parameter values and search for mathematical formulations with which the behavior of the model best matches that of experiments. With automated high-resolution life-cell imaging techniques, large amounts of experimental data are being collected in cell biology labs worldwide. To accelerate the search, the numerical program of this optimization in a large multidimensional space will be parallelized and performed on the high-performance computers with hundreds of computer nodes. As shown in the flowchart in Fig. 1, different sets of model parameters and model equations are used as the input set. The first generation of the sets will be evaluated by their fitness and each computer node runs calculations for each individual in a population by using distributed parallel computing. Characteristic of cell migration, such as cell shape, cytoskeleton architecture, migration speed, etc., will be extracted from the simulations and used to compare with the metrics calculated from the experimental observations. Following that, if the convergence criteria is met, the search will stop, and if not, the calculation will be continued by producing the second generation using evolutionary methods; for example, crossover, mutation, etc. The second generation will be evaluated according to their fitness and so on. The proposed model allows one to identify which set of parameters and equations (i.e., assumptions or hypotheses) will match the specific cell migration behavior. The model and the computer simulation program, once developed and validated, will be used in computer-aided design in tissue engineering. For example, it will be used to design suitable biomaterials with optimal mechanical properties, the optimal topology of ECM, and the 3D spatial placement of cells, to facilitate microtissue formation. In tissue-engineering applications, biological and chemical parameters are frequently considered, while the equally important physical/mechanical design variables have often been neglected. For a rational design of tissue engineering, however, all variables influencing cell function and tissue morphogenesis must be considered. This proposed computational model on microtissue formation will enable the integration of chemical, mechanical, and topographical aspects of the problem and can have a powerful impact on the rational design in 3D bioprinting.

a_para_basic.m
Input variables for the model parameters.

a_driver_durotaxis.m
This driver will set the other required parameters for the durotaxis simulation and within the driver we can either choose to run a simulation or plot the previous simulation results. Therefore the driver calls either a_multicellular_system.m or the plot_simulation_results.m. Note, that there are different drivers for different scenarios.

a_multicellular_system.m
This is the main function of the program, it will call the initialization subroutine and starts the time integration loop, follows the detailed flowchart in Figure 3, and calls for the solvers. It will end the simulation whenever the simulation time reaches.

calc_delta_t.m
This subroutine calculates the time step dt .
The algorithm adopts an adjustable time step to make sure that the time integration is stable. It estimates the maximum possible time step for numerical integration based on the following two criteria: 93 1. dt < 1/10 * min{time scales of decay for all time-dependent variables} 2. dt < (1/4 * (smallest element size) ^2/min{Diffusion constants} After the estimation if the input time step is bigger than the estimated time step. It assigns the estimated time step to the time step variable dt to make sure that we have an stable time integration.

a_multicellular_initialization.m
This subroutine initializes the matrix/tensors/vectors we are using for the finite element simulation. It will call following subroutines:

mesh_a_cell_m.m
This subroutine is for discretizing the 2D cell domain based on the initial cell shape and the element type. Based on the para.initial_cell_shape input, we can switch between several cell shapes and element type. This algorithm also returns the edge segments of the domain.

mesh_find_and_sort_cell_edge.m
The input for this algorithm is the nodal coordinates and the connectivity matrix from the mesh_a_cell_m.m subroutine. And the output is the edge segments.

mesh_brand_new_remeshing.m
94 this subroutine is for the automatic mesh generation. It creates initial distribution in bounding box, add more points near the edge, add some random noise to the coordinates in case the algorithm does not converge well, remove points outside the region, applies Delaunay triangulation to create elements.

mesh_remove_narrow_membrane_tub.m
remove narrow and tube elements. This happens when cells become too narrow at some zones.

mesh_fix_delaunay_mesh.m
this algorithm fixes the meshes after the Delaunay triangulation. Several problems may occur after the default Delaunay triangulation: multiple loops, interior nodes come to edge, and also check if multiple loops exist due to Delaunay triangulation, keep the biggest one.

mesh_clear_singular_node.m
This algorithm removes the singular nodes. In a regular good topology, simple topology, each edge node has two edge segments. Singular nodes have 4 edge segments; we use this rule to find singular nodes

mesh_find_cell_edge_only.m
This subroutine assumes the element nodal numbering is counter clockwise.

mesh_find_edge_normal.m
Finds the normal to the boundary/edges of cell on the nodes.

mesh_Laplacian_smoothing.m
Without changing the connectivity matrix of mesh, optimize the nodal positions to obtain high quality mesh.

apply_zero_FA_condition.m
This subroutine applies adhesive island boundary condition by setting the focal adhesion of those region of the substrate equal to zero.

update_microtubule.m
This subroutine updates the vector associated with microtubule.

update_Stress_Fiber.m
This subroutine updates the stress fiber based on the model.

solver_reaction_diffusion.m
This subroutine solves the system of equations resulted from the finite element model of the reaction diffusion equation. It also has option to switch for different type of the elements.
solver_elasticity.m the displacement. It's include assembly of the mass and stiffness matrix and solving the system of equations.

solver_rate_equation.m
This subroutine is for solving the rate equation. For example, the focal adhesion rate equation.

solver_FEA_K_M.m
This solver pre-calculates the mass (M) and stiffness (K) matrix to avoid calculating that in each time steps.

update_cc_adhesion_bond.m
updates the cell-cell adhesion bond in each iteration.

calc_cell_retraction.m
calculate retraction of the cell using the Eq. (12) get_cell_neighbors.m update the neighbors cell list based on a cut off distance.

calc_P_source.m
97 this subroutine is for calculating the source term of the reaction diffusion equation for the protrusion.

update_variables_after_growth.m
This subroutine updates the variables after the growth happened.

update_variables_after_remeshing_m.m
This subroutine transfers the field variables to the new sets of elements after remeshing.

calc_cell_growth_node_based.m
this subroutine updates the protrusion displacement of the cells based on Eq. (11)

Appendix B (Variables)
The parameter values for the simulations will be chosen from the available experimental data in the literature, the rest of the parameters will be chosen in a fashion to obtain similar results to experimental studies in the literature. Following are the parameter values used for the simulations in this study unless specifically mentioned.