Inverse tissue mechanics of cell monolayer expansion
The body of multicellular organisms must be properly shaped in order to exert its functions , and this proper formation is based on the orchestration of cellular behaviors , such as cell division , differentiation , migration , and other ., One of the key processes in morphogenesis is the coordinated change in cell shapes and positions ., The coordination depends on cell-generated mechanical forces that introduce stress , which induces multicellular deformation and flow 1 ., Therefore , the research on how the molecular components responsible for force generation and propagation , such as motor proteins and cell-cell adhesion molecules , are regulated in space and time during the morphogenesis has attracted a lot of attention recently 2 ., In parallel , remarkable progress has been made in the development of the technologies allowing the measurements of the generated forces and stress in the living tissues 3 , which represents a crucial step towards linking the underlying molecular activities with the morphogenesis ., Epithelial tissues represent important model systems for the understanding of the force dynamics during morphogenesis , because their two-dimensional sheet structure facilitates the observation of the processes that occur in these tissues and analysis ., In particular , many valuable insights have been obtained using the cultured cell monolayer , i . e . , one-cell-thick sheet of tightly-connected epithelial cells 4–6 ., The cells belonging to a monolayer collectively migrate in order to fill a cell-free surface , which replicates in vivo tissue remodeling , such as wound repair , which occurs during regeneration , and epiboly , during embryonic development ., When migrating , the cells exert forces on the underlying substrate to propel themselves forward , and in the unicellular motion , this force , known as the cell traction force , can be visualized by the displacement of fluorescent beads embedded into the substrate 7 ., The simple flat-sheet structure of the monolayer allows us to apply the same technique to observe a spatio-temporal profile of the cell traction force in a wide field of view 8 , and to determine where and how the force and stress are generated 9 , 10 ., In order to achieve the quantitative understandings of the resultant tissue morphogenesis , however , we need to elucidate the other mechanical factors as well , i . e . , the mechanical properties that describe the relation between the deformation and forces ., Although several pioneering works exist 11–14 , our access to the mechanical properties is still limited ., The characterization of these properties often requires exogenous manipulation of the tissue to induce deformation , but the procedure itself perturbs cell physiology and interferes with the tissue morphogenesis ., Here , force measurement in a non-invasive manner gives a way to bypass this issue , and we can infer mechanical properties by associating spontaneous tissue deformation with the observed force dynamics ., In this study , we propose a reverse-engineering method to identify the mechanical properties , which is based on the combination of tissue mechanics modeling and statistical machine learning ., Our strategy is to represent a cell monolayer as a continuum-mechanical system 15 , and to use the passive and simultaneous observations of the deformation and traction force in order to compute the maximum likelihood estimate of the mechanical parameters ., We formulated the inference as an inverse of the forward processes in which the mechanical properties and reaction force to the traction cause the tissue deformation ., Our inference algorithm is based on the sequential updates of estimates; using the current model state and parameters , the mechanical model predicts the traction force field , and then the error feedback based on the observation is used to update the model state and parameters ., Here , we applied our method to a cultured monolayer system to infer the elastic moduli from the collected tissue deformation and traction force data ., To characterize the tissue deformation , we used velocity field of tissue motion , hereafter called tissue flow field ., MDCK cells ( strain II ) were maintained in minimal essential medium ( MEM; Invitrogen ) supplemented with 10% fetal bovine serum ( FBS; Equitech-Bio ) , GlutaMAX ( Invitrogen ) , and 1 mM sodium pyruvate , in a 5% CO2 humidified incubator at 37 C° ., According to a previously published protocol 9 , 48 h before the image acquisition , 3 μl drop of dense cell suspension ( 8 × 106 cells/ml ) was added to each dish containing the gel and 3 ml medium ., Afterward , 3 h before the image acquisition , the medium was replaced by 3 ml CO2-independent medium supplemented with 10% FBS and GlutaMAX ., For the myosin II inhibition , we added blebbistatin ( Sigma Aldrich ) at a final concentration of 25 μM following the replacement of the medium ., Polyacrylamide gel substrates were prepared according to the previously published protocols 8 , 9 ., Briefly , the gel solution was prepared with 3% acrylamide , 0 . 25% bisacrylamide , 0 . 8% ammonium persulfate , 0 . 08% TEMED ( Bio-Rad products ) , and 0 . 01% red fluorescent carboxylate-modified beads ( 0 . 5μm diameter , Invitrogen ) ., 20 μl of this mixture was added to each dish and the samples were covered with glass cover slips with 18 mm diameter ( Matsunami ) ., After the polymerization , the surface was coated with type I collagen ( Purecol , Advanced BioMatrix ) using 4 μM sulphosuccinimidyl-6- ( 4-azido-2-nitrophenylamino ) hexanoate ( Sulfo-SANPAH; Pierce ) ., Young’s modulus of the gel was characterized by the conventional method using the Hertz equation 16 , obtaining E = 2500±600 Pa ., The Fourier-transform traction microscopy 8 was used to estimate traction force fields from bead displacement fields ., Confocal imaging was conducted at 48 h after the seeding of the cells ., We used FV10i-LIV ( Olympus ) to simultaneously acquire phase contrast images of the cells and fluorescent images of the beads ., The trial period lasted for 6-10 h and the sampling rate was one frame per 5 min ., After each trial , we removed the cells by the trypsinization and imaged the strain-free pattern of the fluorescent beads ., To increase the field of view , we stitched tiled images by the Grid/Collection stitching plugin in Fiji 17 ., Following this , the images at different time points were aligned to match the bead configurations in a cell-free region 18 ., To obtain velocity fields in the phase contrast image and bead displacement fields , we adopted an advanced optical flow technique , which tracks changes between two images by matching the patterns of intensity and its gradient 19 ., S1 Movie shows a representative result of the image analysis ., For the tissue flow , we used images from subsequent time points , while for the bead displacement , we compared the stress-free image with each fluorescent image ., The image resolution was 0 . 61μm/pixel and the grid spacing of the vector fields was 14 . 7 μm ., Finally , the flow and force fields were down-sampled in space and time into Δx = 29 . 4μm and Δt = 10 min ., We adopted a continuum modeling of the monolayer mechanics , in which the deformation of the tissue , or strain , determines the stress 20 , 21 ., The cell monolayer was represented as a two-dimensional sheet , and therefore , we represented the stress as a symmetric matrix, σ ( x , y , t ) ≡ ( σ x x σ x y σ x y σ y y ) = π I + σ ˜ ., ( 1 ), In the second line , we applied the deviatoric decomposition where the stress tensor is given as the summation of isotropic ( the first term ) and distortional ( the second term ) components ., The strain tensor was also represented by a two-by-two matrix as, ϵ i j ≡ 1 2 ( ∂ j u i + ∂ i u j ) , ( i , j = x , y ) , ( 2 ), where ( ux , uy ) represents the displacement vector of the tissue from the stress-free state ., In a linearly elastic material , the relationship between the stress and strain tensors becomes simply linear , meaning that the stress accumulates in response to the strain ., However , the stress-free state of the living tissue can vary in time due to cell growth and death ., Therefore , we adopted an alternative formulation using the strain rate tensor, e ( x , y , t ) ≡ ϵ ˙ → e i j = 1 2 ( ∂ j v i + ∂ i v j ) → e = 1 2 ( ∇ · v ) I + e ˜ , ( 3 ), where v = ( vx , vy ) is the flow velocity vector in the tissue ., In the last line , we applied the deviatoric decomposition ., Although previous works suggested that anisotropic cell division can contribute the tissue mechanics 22 , 23 , we modeled the cell growth simply as isotropic and homogeneous expansion with the rate Dg , which is partially supported by a previous report that cell division in the monolayer shows no particular orientation 24 ., Since the observed total expansion of the tissue is the summation of the growth and the deformation-originated expansion , i . e . , Dtotal = Dg + Dmaterial , the subtraction Dtotal − Dg should appear in the stress-strain relation 15 ., Taken together , our elastic model was written as, π ˙ ( x , y , t ) = K ( ∇ · v ( x , y , t ) - D g ) + ξ ( x , y , t ) σ ˜ ˙ x x ( x , y , t ) = 2 G e ˜ x x ( x , y , t ) + ξ x x ( x , y , t ) σ ˜ ˙ x y ( x , y , t ) = 2 G e ˜ x y ( x , y , t ) + ξ x y ( x , y , t ) , ( 4 ), where K and G are the in-plane bulk and shear elastic moduli , respectively ( in S1 Text , we derived the relation of the in-plane moduli to the conventional three dimensional moduli ) ., ξs are the stochastic terms representing random variables with Gaussian distribution that is not space or time-dependent ., Dg is associated with the cell division interval tdiv as Dg = ln 2/tdiv , and we adopted tdiv = 1 division per day ., We found that essentially the same results are obtained by increasing or decreasing the rate of growth rate two times ., On the other hand , the tissue stress tensor and traction force vector were related through the force balance equation 25, - T x ( x , y , t ) = ∂ π ∂ x + ∂ σ ˜ x x ∂ x + ∂ σ ˜ x y ∂ y + η x - T y ( x , y , t ) = ∂ π ∂ y - ∂ σ ˜ x x ∂ y + ∂ σ ˜ x y ∂ x + η y , ( 5 ), where ηs are noises in the force quantification assumed to be normally-distributed , and we call them the observation noises ., Here , we briefly describe the inference algorithm ( the details of derivation is given in S1 Text ) ., Let Y and Λ represent the collected spatio-temporal fields of traction force and tissue flow , respectively , and X represent the stress tensor field that was discretized in space and time according to Y and Λ ., Additionally , let θ represent the model parameters ., Then , our aim is to find such θ ^ that maximizes the log-likelihood:, ln L ≡ ln p ( Y | Λ , θ ) = ln ∫ p ( Y | X , θ ) p ( X | Λ , θ ) d X ., ( 6 ), Note that p ( X|Λ , θ ) and p ( Y|X , θ ) are corresponding to the stress evolution and force balance equations , i . e . , Eqs 4 and 5 , respectively ., Since the integration w . r . t . X is analytically intractable , we adopted the expectation-maximization ( EM ) algorithm , which maximizes the lower-bound of the log-likelihood by executing the following E and M steps alternately 26 ., ( E-step ), Estimate the stress fields by computing p ( X|Y , Λ , θ* ) with the Rauch-Tung-Striebel smoother 27 , where θ* is a tentative estimate of the parameters ., ( M-step ), Compute the expected complete-data log-likelihood:, Q ( θ ) = E p ( X | Y , Λ , θ * )  ln p ( X , Y | Λ , θ )  , ( 7 ), and update the parameters through maximizing Q ( θ ) ., Repeating the E and M steps , which offers monotonic increase and convergence of the likelihood ., After the convergence , we obtain the maximum likelihood estimate of the parameters θ ^ ., In order to collect the data on tissue deformation and force , we adopted a model system , Madin-Darby canine kidney ( MDCK ) epithelial cell monolayer , and analyzed it using the colony expansion assay 9 ., We performed phase contrast imaging to measure the flow of cells in the monolayer and , simultaneously , traction force microscopy , in order to visualize the generated force using fluorescent beads embedded into the soft substrate ., Our inference algorithm used a mechanical model of the cell monolayer , i . e . , a spatio-temporal model of mechanical stress within the tissue ., Our mechanical model of the cell monolayer , represented by Eq 4 ( see Materials and methods section ) , included two biophysical factors that are essential for the colony expansion: tissue elasticity and cell growth ., Elasticity is a basic property of a material , which resists the influence of an external force and shows a recoverable deformation ., According to the previous studies 9 , 14 , we assumed the linear elasticity where the deformation is proportional to the force ., Additionally , cellular growth supplies new cells into the tissue , and thereby promotes tissue expansion 28 ., Our mechanical model also included stochasticity , which represents other mechanical processes such as viscosity , plasticity , active contractile force , and others ., This model , represented by Eq 4 , had two mechanical parameters describing the elastic properties of the monolayer , the in-plane bulk modulus K and shear modulus G . The values of the bulk and shear moduli represent the resistance against area-changing and area-preserving deformation in the monolayer , respectively ., Additionally , the elastic model contains the variance parameters for strength of the stochastic effects in the stress dynamics ., Our inference algorithm , using the movie data showing tissue flow and traction force , estimated the values of these parameters ( Fig 1A and 1B and “Materials and methods” section ) ., We found that the flow speed of the tissue at the periphery was approximately 10-30 μm/h ( Fig 2A ) , the strength of the traction force was distributed around 10-100 Pa ( Fig 2B ) , and both flow speed and force strength decreased monotonically along the distance from the edge ( Fig 2A and 2B ) ., These results are consistent with the previous observations 4 , 8 , 29 ., We computed maximum likelihood estimates of the parameters in our elastic model from the collected data ., For this estimation , the model state , i . e . , the tissue stress field , was corrected by the current traction force data at each time point in the movie sequence; following this , the model was numerically simulated , using the tissue flow data , in order to predict the traction force field at the following time point ., As a result , even though the model inference was based on the one-step prediction ( Δt = 10 min ) , we found that the estimated model can provide a long-term forecast ( >1 h ) without corrected by the traction force data ., To quantitatively demonstrate this result , we divided each movie data on tissue flow and traction force , into two parts in time: The earlier , training data and the following test data ., Using the training data , the inference algorithm was used to estimate the model parameters and the stress field in the monolayer , and then this model was examined in terms of the forecast accuracy for future force fields by using the test data ., As a quantitative measure , we employed the correlation between the forecasted and observed force vector fields:, R = ⟨ T forecast · T data ⟩ ⟨ | T | forecast · | T | data ⟩ , ( 8 ), where 〈⋅〉 represents an average over all spatial grid points ., The correlation plotted against time is represented in Fig 3A ., As shown , the forecast provided by the elastic model was highly correlated even 3 h after the initiation of the test part ., For comparison , we adopted a null-hypothetical , zero-elasticity model , where K = G = 0 ( Fig 3A ) ., In Fig 3B , the correlation in both models at the last time-point in the test part , the long-term forecast accuracy , is shown ., These results demonstrate that our data-driven elastic model showed better forecast which is clearer especially in longer time forecast ., We also computed the difference of correlation between the models in a sample-wise manner ( Fig 3C ) , and confirmed statistically-significant superiority of the elastic model compared with the zero-elasticity model ( p < 10−4 , Wilcoxon signed-rank test ) ., Representative forecasted and observed traction forces at the last time point are shown in Fig 3D ., Therefore , despite of its simplicity , our data-driven elastic model captured the stress evolution in the tissue expansion by the estimated bulk ( K ) and shear ( G ) moduli ., Next , we examined if the elastic moduli are different in the tissues treated with blebbistatin , a myosin inhibitor ., Previous studies showed that the inhibition of the molecular motors considerably reduces the traction force strength , which is expected ., However , this inhibition does not slow down the tissue expansion rate 6 , indicating the alterations in tissue mechanical properties ., By comparing the elastic moduli estimated from a different dataset from blebbistatin-treated tissues with those estimated under the standard conditions ( Fig 4A , 4B and 4C ) , we obtained the results consistent with those obtained in previous experiments ., We observed a significant decrease in the elastic moduli associated with the treatment ( Fig 4D and 4E ) ., In the standard experimental setting , both moduli were within the order of the magnitude of ∼ 103Pa · μm , while myosin inhibition induced several-fold reduction in the elastic modulus values , i . e . , softening ( p < 0 . 01 , U-test ) ., Note that the estimated moduli are not guaranteed to have positive values because unmodeled monolayer mechanics , such as viscosity and anisotropic tissue growth , might affect the stress dynamics ., In fact , the estimated moduli from myosin-inhibited tissues have frequently shown negative values , suggesting that the elasticity effect was no longer dominant over unmodeled effects due to the softening ., Finally , we assessed the spatial distribution of the elastic moduli ., When considering difference in the flow speed and force strength dependent on the distance from the tissue edge ( Fig 2 ) , we expected the mechanical properties that would correlate with this distance ., However , as shown in Fig 5 , we found that the estimated moduli are homogenous along with the distance ., The framework we presented in this study would benefit from the advances in force measurement ., For example , although the traction force microscopy is applicable only to in vitro tissues , in vivo measurement techniques are being actively developed 38 , 39 ., We can apply our reverse-engineering method to the in vivo measuring by modifying the model of the observation process , Eq 5 in our case ., Additionally , another interesting direction would be to use a more elaborate model of tissue mechanics , in particular , by directly including cellular processes such as cell division 40 , 41 ., We hope that , with the advancements in the technology of force/stress measurement , our method may assist further understanding of the mechanics underling tissue development and maintenance .
Living tissues undergo deformation during morphogenesis ., In this process , cells generate mechanical forces that drive the coordinated cell motion and shape changes ., Recent advances in experimental and theoretical techniques have enabled in situ measurement of the mechanical forces , but the characterization of mechanical properties that determine how these forces quantitatively affect tissue deformation remains challenging , and this represents a major obstacle for the complete understanding of morphogenesis ., Here , we proposed a non-invasive reverse-engineering approach for the estimation of the mechanical properties , by combining tissue mechanics modeling and statistical machine learning ., Our strategy is to model the tissue as a continuum mechanical system and to use passive observations of spontaneous tissue deformation and force fields to statistically estimate the model parameters ., This method was applied to the analysis of the collective migration of Madin-Darby canine kidney cells , and the tissue flow and force were simultaneously observed by the phase contrast imaging and traction force microscopy ., We found that our monolayer elastic model , whose elastic moduli were reverse-engineered , enabled a long-term forecast of the traction force fields when given the tissue flow fields , indicating that the elasticity contributes to the evolution of the tissue stress ., Furthermore , we investigated the tissues in which myosin was inhibited by blebbistatin treatment , and observed a several-fold reduction in the elastic moduli ., The obtained results validate our framework , which paves the way to the estimation of mechanical properties of living tissues during morphogenesis .
In order to shape the body of a multicellular organism , cells generate mechanical forces and undergo deformation ., Although these forces are being increasingly determined , quantitative characterization of the relation between the deformation and forces at the tissue level remains challenging ., To estimate these properties , we developed a reverse-engineering method by combining tissue mechanics modeling and statistical machine learning , and then tested this method on a common model system , the expansion of cultured cell monolayer ., This statistically sound framework uses the passive observations of spontaneous deformation and force dynamics in tissues , and enables us to elucidate unperturbed mechanical processes underlying morphogenesis .
tissue mechanics, fluorescence imaging, mechanical properties, classical mechanics, cell cycle and cell division, cell processes, biomechanics, developmental biology, molecular motors, actin motors, materials science, damage mechanics, morphogenesis, motor proteins, research and analysis methods, contractile proteins, imaging techniques, proteins, deformation, biophysics, physics, biochemistry, cytoskeletal proteins, cell biology, myosins, biology and life sciences, physical sciences, material properties