Mechanism underlying dynamic scaling properties observed in the contour of spreading epithelial monolayer
MMechanism underlying dynamic scaling properties observed in the contour ofspreading epithelial monolayer
Toshiki Oguma, Hisako Takigawa-Imamura, and Takashi Miura
Department of Anatomy and Cell Biology, Graduate School of Medical Sciences, Kyushu University, Japan (Dated: July 6, 2020)We found evidence of dynamic scaling in the spreading of MDCK monolayer, which can be char-acterized by the Hurst exponent α = 0 .
86 and the growth exponent β = 0 .
73, and theoreticallyand experimentally clarified the mechanism that governs the contour shape dynamics. During thespreading of the monolayer, it is known that so- called ”leader cells” generate the driving force andlead the other cells. Our time-lapse observations of cell behavior showed that these leader cells ap-peared at the early stage of the spreading, and formed the monolayer protrusion. Informed by theseobservations, we developed a simple mathematical model that included differences in cell motility,cell-cell adhesion, and random cell movement. The model reproduced the quantitative characteris-tics obtained from the experiment, such as the spreading speed, the distribution of the increment,and the dynamic scaling law. Analysis of the model equation revealed that the model could repro-duce the different scaling law from α = 0 . , β = 0 .
25 to α = 0 . , β = 0 .
75, and the exponents α, β were determined by the two indices: ρt and c . Based on the analytical result, parameter estimationfrom the experimental results was achieved. The monolayer on the collagen-coated dishes showeda different scaling law α = 0 . , β = 0 .
68, suggesting that cell motility increased by 9 folds. Thisresult was consistent with the assay of the single-cell motility. Our study demonstrated that thedynamics of the contour of the monolayer were explained by the simple model, and proposed a newmechanism that exhibits the dynamic scaling property.
I. INTRODUCTION
The shape of mammalian cell colonies varies, depend-ing on the cell type and its environment. In the fieldof oncology, there is known to be a correlation betweenshape and malignancy of cancer, and suitable strategiesfor treatment can be inferred from analyzing the shapeof cancer cell colonies [1]. To quantify the shape of thesecolonies, fractal analysis is often used. It has been re-ported that a high fractal dimension D reflects a hetero-geneous contour shape, and that fractal dimension andcancer malignancy are likewise correlated [2]. Cancercells show higher proliferation rates, higher motilities,and weaker cell-cell adhesion than benign cells [1]. Thesedifferences are thought to affect collective cell behavior,and contribute to the roughness of contour shapes in can-cer.Cell movements are promoted by chemical or physi-cal cues such as signal molecules and mechanical forces,in response to which cells modulate their downstreamcytoskeleton by altering the subcellular localization ofsmall G proteins such as RhoA, Rac1, and Cdc42 [3–5].In collective cell migration, highly motile cells are oftenobserved at the edge of the epithelial monolayer in vitro .Known as leader cells, these are characterized by havinghigh Rac1 activity, thick actin filaments, and large cellbodies with spreading lamellipodia. Leader cells can gen-erate a driving force to spread the epithelial tissue, andare often observed at the tips of monolayer protrusions[6–8]. However, the relationship between leader cells andthe overall contour shape remains controversial. Markand colleagues suggested that sharp curvature promotedthe appearance of the leader cells [9]. On the other hand,several studies have suggested that leader cells are deter- mined by different mechanisms, such as mechanical force[10] and the Dll4-Notch1 pathway [11], then form mono-layer protrusions.It is known that many curves in nature, such asthe earth’s surface in cross-section and the interface ofclouds, form self-affine fractal structures [12]. The frac-tal dimension D is used to quantify the roughness of thestructure. It can be measured using the box-countingmethod [13, 14]. For example, the simplest self-affinefractal structure is a Brownian curve, which is the tra-jectory of a particle undergoing Brownian motion plot-ted over time, and its fractal dimension is D = 1 .
5. Ina self-affine fractal, a similar curve can be obtained byexpanding a smaller part of the curve.If we let a and b , the magnification rates for scaling,correspond to coordinate axes in (1 + 1) dimension, thepower exponent α , which satisfies a = b α , is called aHurst exponent. Self-affine fractal structure is character-ized by the Hurst exponent α , and it is known that α and the fractal dimension D satisfy α + D = 2 [13, 15].For a growing interface, the local roughness w ( l, t ) isdefined as the standard deviation of the height of the con-tour, within the closed range l . Such a system is said toexhibit dynamic scaling when the following relationshipsare satisfied [14]: w ( l, t ) ∼ (cid:34) l α for l (cid:28) l ∗ t β for l (cid:29) l ∗ , (1)where the length l ∗ increases as time evolves. The formerequation confirms that the interface forms a self-affinefractal structure, and the latter means that the hetero-geneity of the whole interface is scaled by the time. Theexponent α is the Hurst exponent of the interface, and β is called the growth exponent. a r X i v : . [ q - b i o . CB ] J u l Dynamic scaling has been identified in a wide range ofphenomena, including bacteria colonies on agar [16, 17],propagation of slow combustion of paper [18], and grow-ing interfaces in a turbulent liquid crystal [19]. Incollective migration of vertebrate cells, the contoursof cell colonies of HeLa and Vero cell lines were alsofound to exhibit dynamic scaling, with α = 0 .
50 and β = 0 .
32 [20, 21]. Some partial differential equationsthat show the properties of dynamic scaling are well-known, including the Edwards-Wilkinson (EW) equa-tion, the Kardar-Parisi-Zhang (KPZ) equation, and theKuramoto-Sivashinsky (KS) equation. The EW equationis a stochastic partial differential equation with diffusionand spatiotemporally independent noise terms [22], theKPZ equation is similar to the EW equation but includesa non-linear term [23], and the Kuramoto-Sivashinskyequation is a partial differential equation with the dif-fusion, non-linear, and spatially fourth derivative terms[24]. Intriguingly, both the KPZ and KS equations areknown to have the same Hurst and growth exponents, α = 1 / β = 1 /
3. The sharing of similar scalingexponents between different phenomena and solutions isa property known as universality. In particular, the uni-versality characterized by α = 1 / β = 1 / α and β are determined by the index c , which reflects therelative intensity of the random movement on differencesof cell motility. II. MATERIALS AND METHODSA. MDCK cell culture
MDCK II cells were cultured using Eagle’s minimalessential medium (MEM; Nacalai tesque) containing 10% fetal bovine serum (FBS) and 100 U/ml penicillin-streptomycin (Nacalai tesque), and maintained in a 5 %CO controlled atmosphere at 37 ◦ C. B. Fabrication of PDMS sheets
Polydimethylsiloxane (PDMS) sheets, used to create acell-free regions, were prepared in the following manner:A well-mixed PDMS (Sylgard 184, Toray) precursor so-lution, with a 10:1 ratio of prepolymer to curing agent,was poured onto a flat polystyrene plate to a thicknessof 1 mm, then cured in a drying oven at 80 ◦ C for 2 hr.After curing, 8-mm PDMS disks with a 3- or 4-mm holeswere created using biopsy punches (Maruho).
C. Cell patterning
The PDMS sheets were placed on the surface of a27 mm glass bottom dish (IWAKI), and the holes inthe sheets were filled with MEM (without air bubbles),along with 2 ml of MDCK cells suspended in medium(2 . × cells/ml). Samples were incubated at 37 ◦ Cand 5 % CO for 48-72 hr, until cells reached conflu-ence within the holes, then the PDMS sheets were gen-tly removed and washed with MEM twice. The mediumwas changed with fresh medium containing CellTrackerGreen CMFDA Dye (5 µM ; Invitrogen) and Hoechst33342(1 µg/ml ; Dojindo). After 4 hr incubation, thetime-lapse observation was performed. The collagen-coated dishes were prepared with 27 mm glass bottomdishes and type I-c collagen solutions (100 µg/ml ; NittaGelatine). D. Time-lapse microscopy
The time-lapse observations were performed using aNikon A1 confocal microscope with a 10 × or 20 × ob-jective lens. The cells were maintained in a 5 % CO controlled atmosphere at 37 ◦ C. Images were acquiredevery 20 min (for Fig. 3) or 30 min (for Fig. 2, 6, S5)until 18 hr had passed after the PDMS sheet removal.
E. Measurement and quantification of cellspreading
The green fluorescent microscopy images were numer-ically converted and analyzed quantitatively using Im-ageJ (NIH), Python, and Mathematica (Wolfram). Theimages were binarized to capture the shape of the mono-layer, with the threshold values for binarization deter-mined and set manually. We defined the centroid as thecenter of gravity of the cell monolayer in the initial im-age. The distance D ( θ ) from the centroid to the contourwas measured along 2000 directions, with constant inter-vals. Then, the mean front distance was calculated as (cid:104) D ( θ ) (cid:105) ( t ) = (cid:80) θ D ( θ ) / w ( l, t ) and l within the range of l ∈ [2∆ x, x ].Here, ∆ x was defined by 2 π (cid:104) D ( θ ) (cid:105) / w ( l max , t ) and t . The values shownin Figs. 2(b,c,e), 6(a), and S5(c) are averages of the ex-perimentally obtained values. F. Cell tracking
We labeled cell nuclei with Hoechst33342, and imageswere obtained from 1 hr to 18 hr after PDMS removal andanalyzed with ImageJ (NIH), using the TrackMate plu-gin [26] to manually track the centers of cell nuclei. Weconsidered cells whose nuclei were located within 30 µm of the contour to be at the edge of the monolayer. Thecontours were determined from the brightfield images,and manually traced with segmented lines. G. Cell staining
MDCK cells were washed with phosphate-bufferedsaline (PBS), fixed with 4% paraformaldehyde (PFA)for 10 min, and permeabilized in 0.1 % Triton X-100/PBS for 10 min. After washing, cells were stainedwith Alexa Fluor 555 Phalloidin (1:40; Invitrogen) andHoechst33342 (1:2000; Dojindo). Cells were incubatedfor 30 min at room temperature before observation. Flu-orescent and phase-contrast images were acquired usinga BZ-X810 fluorescence microscope (Keyence) with 20 × objective lens. H. Numerical simulation
The numerical simulations were performed usingMathematica and Julia, and we used the explicit Eulerscheme for calculating the time evolution. The code usedfor this study is available from the corresponding authorsupon request.
I. Single cell tracking assay
MDCK cells were seeded onto the 27 mm glass bottomdish (2 . × cells per ml). The cells were incubated at37 ◦ C and 5 % CO for 24 hr until the cells adhered to thebottom of the dish. Hoechst33342 (1:2000) was added tothe dish and incubated for 30 min. Images were acquiredusing a Nikon A1 confocal microscope at 5 min intervalsover a span of 2 hr. The centers of cell nuclei were man-ually tracked using the TrackMate plugin [26] in ImageJ,and statistical analysis was performed in Mathematica,using Student’s t-tests to compare experimental groups. III. RESULTSA. Epithelial monolayer spreading
To investigate how the epithelial monolayer spread,MDCK cells were cultured in the closed circular areaconfined within the PDMS sheet. Time-lapse observa-tions were performed after removal of the PDMS sheetboundary (Fig. 1(a), Supplemental Movie 1). (a) PDMS sheet8 mm4 mm MDCK cells(b) Contour ExtractionRaw data PDMS sheet removal (0 hr) 18 hr
FIG. 1. Overview of the experimental method and mea-surement procedures. (a) Observation of MDCK migration.MDCK cells were seeded in the region confined by the PDMSsheet. After cells reached confluence, the PDMS sheet was re-moved and time-lapse observations were performed. (b) Mea-surement of the contour shape. The MDCK monolayer wasvisualized by Cell Tracker Green, and the microscopy imageswere transformed into binarized images. The distances D ( θ )from the center to the edge of the monolayer were measured.The orientation of the points along the contour was definedas θ (0 < θ < π ). The angle φ = l/ (cid:104) D (cid:105) ( t ) was defined fromthe unit of measurement l . Scalebar = 1000 µm . While the contour of the epithelial monolayer wasinitially smooth and round, our observations revealedthat it became rougher and more uneven over time, ascells migrated to the cell-free area. The epithelial cellskept contact each other (Supplemental Movie 1). Toquantify this, we measured the distance from the cen-ter to the edge of the monolayer, D ( θ ), where θ indi-cates the orientation of measurement points along thecontour (0 < θ < π ) (Fig. 1(b)). We chose 2000 pointsfor measuring D ( θ ), to accurately reproduce the contourboundary during migration. The initial circumference(after the PDMS removal) was 1 . × µm . We foundthat the average (cid:104) D ( θ ) (cid:105) increased linearly, at a rate of11 . µm/hr , and the standard deviation of D ( θ ) also in-creased (Fig. 2(b) and 2(c)). These results suggestedthat the epithelial monolayer spread at a constant speedwith increasing heterogeneity of the contour. (a) (e)(b) Slope = 0.725Slope = 11.2 (μm/hr) (c)(d) Slope = 0.8584 hr 18 hr4 hr11 hr18 hr FIG. 2. Dynamic scaling law in the spreading of MDCKmonolayer, obtained from experimental observations. (a)Contour of the MDCK monolayer, as it evolved from the ini-tial smooth shape to the rough shape. (b) Time evolution ofthe average values of D ( θ ). Dots indicate the measured val-ues, with the fitted line shown as dashed (n=6). (c) Time evo-lution of w ( l max , t ), the standard deviation of D ( θ ) (n=6). (d)Measurement of Hurst exponent. The plot shows the experi-mental results at different time points, and the black dashedline is the fitted line for 2∆ x ≤ l ≤ x ( µm ) at 18 hr. Fromthis, the Hurst exponent was calculated as: α = 0 . w ( l max , t ) and t . Dots represent the measuredvalues, and the black dashed line is the fitted line (n=6). Thegrowth exponent was calculated as: β = 0 . We characterized the local roughness w ( l, t ) of the cir-cular geometry as follows: w ( l, t ) = (cid:28)(cid:114)(cid:68) [ D ( θ, t ) − (cid:104) D (cid:105) φ ] (cid:69) φ (cid:29) (2) φ = l (cid:104) D (cid:105) , where t (hr) is the time after PDMS removal, and (cid:104)· · · (cid:105) φ denotes the average value in the area φ . The local rough-ness w ( l, t ) is the standard deviation of the distance D ( θ )in the range of φ .To examine whether the power law in (1) holds forthe cell spreading, log-log plots were used to assess therelationship between local roughness w ( l, t ) and l at dif-ferent time points (Fig. 2(d)). Within the range of small l , ln w ( l, t ) linearly increased along with increasing ln l ,indicating that the power law w ( l, t ) ∼ l α holds, whichsuggests the contour of the epithelial monolayer is a self-affine fractal structure. We also found the range of l that showed a linear relationship to expand over time, whichis consistent with the well-known observation that l ∗ in-creases with time under the dynamic scaling law [14, 25].We determined the Hurst exponent to be α = 0 . t = 18. For large l , the value of w ( l, t ) approaches w ( l max , t ), where l max = 2 π (cid:104) D (cid:105) . As shown in Fig. 2(e),ln w ( l max , t ) and ln t show a linear relationship, indicatingthat the power law w ∼ t β holds. The growth exponentwas β = 0 . α = 0 . ± .
013 and β = 0 . ± . B. Cell behavior at the edge
To understand the dynamics of the evolving contourshape as the monolayer spreads, we observed and quan-tified the behavior of cells near the monolayer edge. Cellnuclei were visualized and tracked from t = 1 to t = 18(Fig. 3(a), Supplemental Movie 2). Among cells initiallylocated in the edge region, 58% of these remained nearthe edge the entire observation time, while 42% migratedtowards the monolayer center. These internalizationswere observed to result from the merging of monolayerprotrusions. The decrease in the number of cells aroundthe edge was compensated for by proliferation and theintercalation of internal cells. We reversely tracked cellsnear the edge at t = 18, and found that 77 % of thesewere also near the edge at t = 1. In addition, the kymo-graph of the edge cell migration showed that movementtowards the initially cell-free region first started amongedge cells, then the inner cells followed (Fig. S1). Theseresults suggested that the dynamics of the contour shapewere primarily driven by the movements of cells at ornear the edge.Since leader cells are known to play an important rolein the formation of monolayer protrusion, we next fo-cused on their behavior and dynamics. Leader cells arecharacterized by their large lamellipodia, large cell bod-ies, and high motility. Cells with these characteristicshapes were observed as early as t = 2 (Fig. 3(b)). Asshown in Fig. 3(c), we performed reverse tracking of theleader cells from t = 18 to t = 1. Leader cells were clearlydistinguished by their locations and large cell bodies at t = 18. They were consistently at the edge, formed themonolayer protrusions, kept located at the tip, and hadhigher velocities than other cells (Fig. S2). The fact thatleader cells kept located at the tip of the protrusion re-flected its spontaneous high motility, and suggested thatthe high motilities of the leader cells were maintainedthroughout the observation period. These results indi-cate that leader cells emerged at an early stage of themigration at the edge, and that their properties did notchange. FIG. 3. (Color) Behavior of the cells at the edge of the monolayer. (a) Cell tracking at the monolayer edge. Nuclei werevisualized with Hoechst33342 (left panels, blue circles). Cells that remained near the edge from t = 1 to t = 18 are indicatedby red circles. Those cells that dropped and intercalated from the edge by t = 18 are indicated by green and yellow circles,respectively. Among 45 cells initially located at the edge, 19 cells dropped, while 11 cells divided (right panel). 11 cells newlyintercalated to the edge, therefore, 77 % of the edge cells at t = 18 are originated from the cells at the edge at t = 1. Scalebar = 100 µm . (b) Snapshot of the edge at t = 2. Large cells that stretch their lamellipodia were observed (indicated withyellow arrows). Upper panel is a phase-contrast image, and the lower panel is a fluorescent image; the nuclei and F-actin arerepresented in blue and red, respectively. (c) Behavior of leader cells during collective cell migration. Circles indicate leadercells, which are distinguished by their locations and large sizes. Scale bar = 200 µm . C. Mathematical model and numerical simulation
Informed by the experimental results, we modeled thedynamics of the contour of the MDCK monolayer. Weassumed that the contour dynamics arise from the move-ments of cells at the edge, and these cells have differentmotility. The model also assumes, based on a knownproperty of epithelial cells, that cells interact throughintercellular adhesion (Fig. 4(a)). We described thedynamics of cells at the edge by means of temporally-continuous and spatially-discrete differential equations asfollows: ddt r i ( t ) = M i + ( T i − T i + ) + ση i ( t ) r i / | r i | , (3) where r i ( t ) represents the coordinates of the cell i at time t , M i is active, directional movement, and T i describespassive movement due to tension from cell-cell adhesion.The third term on the right-hand side represents randommovement, with the constant σ indicating the intensity ofnoise, while η i ( t ) is spatially and temporally independentGaussian white noise. At the cellular scale, the inertialforce may be ignored.The directional motility term M i was given by M i = (cid:20) v l r i / | r i | (if the cell i is a leader cell) v f r i / | r i | (otherwise) , (4)where v l and v f are the velocity of the leader and followercells, respectively ( v l > v f > p l , and the distribution of leader cellswas randomly determined. (d) (e)(b) (c) Slope = 0.729Slope = 0.849(a) Slope = 11.2 (μm/hr)4 hr 11 hr18 hr FIG. 4. Numerical simulation of the monolayer spreading.(a) Diagram of the mathematical model. Circles representthe cells at the edge, with the dark circles representing leadercells, and the light circles follower cells. The coordinates ofthe cell center are denoted by r i . Cells at the edge have differ-ent types of motility M i , and interact with neighboring cellsthrough cell-cell adhesion T i (b-e) Results from the numericalmodel corresponding to Fig. 2(b) to 2(e). (b) Time evolutionof the average of the distance D ( θ ). (c) Time evolution ofthe standard deviation of D ( θ ) (global roughness). (d) Log-log plot of the local roughness w ( l, t ) against l . (e) Log-logplot of the global roughness w ( l max , t ) against t . Parameters: v l = 94 [ µm h − ], v f = 2 [ µm h − ], p l = 0 . ρ = 6 . h − ],and σ = 3 [ µm h − / ]. Assuming that intercellular tension T i linearly in-creases with the intercellular distance, T i was given by T i = ρ ( r i − r i − ) , (5)where ρ is the tension coefficient.We set the initial shape of the monolayer as a circlewith a radius of 2000 µm . The total number of cells atthe edge was set to N = 2000, aligned with a regularinterval. For simplicity, the cell number change and re-arrangement at the edge were omitted. The numericalresults are shown in Fig. 4 and Supplemental Movie 3.The parameter p l was estimated from the actual distri-bution of cells displaying large lamellipodia (Fig. 3(b)),and v l and v f were set such that the growth speed of theaverage diameter was consistent with the experimentallyobserved values.We found that the model (3) described and capturedthe process of cell sheet expansion well, as it was able toreproduce many of the properties observed in our experi-ments, as shown in Fig. 4(b) to 4(e). The time evolutionof (cid:104) D ( θ ) (cid:105) was linear, with a slope of 11 . µm/h . Thevalues of w ( l max , t ) shown in Fig. 4(c) were similar to those in Fig. 2(c). The log-log plot of w ( l, t ) against l showed a linear relationship for small l (Fig. 4(d)), andthe Hurst exponent was calculated as α = 0 . w ( l max , t ) against t also showed a linearrelationship (Fig. 4(e)), and the growth exponent wascalculated as β = 0 . D ( θ )], Min[ D ( θ )] and the distribution of theincrement for θ were also similar to the experimental re-sults (Fig. S3). Taken together, these results show thatthe model (3) explained the dynamic scaling law seen inthe contour of the MDCK monolayer. D. Analysis of mathematical model
In this section, we show that the Hurst and growthexponent were analytically estimated, then evaluate theeffects of the cell-cell adhesion, the difference of cell motil-ity, and the noise intensity in this system.Assuming that N is sufficiently large, it can be re-garded that the tension affects only the radial direction,thus the model was simplified to a one-dimensional flatmodel with periodic boundary as follows: ddt h x ( t ) = f x + ρ ( h x +1 ( t )+ h x − ( t ) − h x ( t ))+ ση x ( t ) . (6)Here h x ( t ) is the distance from the center to the cell x (for simplicity, set h x (0) = 0), and f x is the motility ofthe cell x that takes v l or v f . The cell with f x = v l wasselected randomly with the ratio of p l . The numericalcalculation of the flat model (6) reproduced the charac-teristic dynamics with the dynamic scaling law generatedby the circular model (3) (Fig. S4).The exponents in the model (6) take different valuesdepending on the parameters. When v l = v f , the model(6) is essentially the EW model, which includes the Gaus-sian white noise and the diffusion term. It is knownthat the EW model shows dynamic scaling α = 0 . β = 0 .
25 [14, 22]. On the other hand, if σ = 0,the model (6) is regarded as the model with temporallyfixed noise and diffusion, which shows the dynamic scal-ing law: α = 0 . β = 0 .
75. Tentatively, we calledthese dynamics the fixed noise model.Considering the discrete Fourier expansions of h x ( t ), f x , and η x ( t ), the Fourier coefficients are denotedby ˆ h k ( t ), ˆ f k , and ˆ η k ( t ), respectively. The randomvariable ˆ f k follows Gaussian distribution denoted as( v l − v f ) (cid:112) p l (1 − p l ) ˆ N (0 ,
1) (See the supplemental textA). The complex Gaussian white noise ˆ η k ( t ) satisfies (cid:82) ts ˆ η k ( t (cid:48) ) dt (cid:48) = ˆ N (0 , t − s ), where ˆ N (0 ,
1) is complex ran-dom variable that follows N (0 , ) + i N (0 , ) (See thesupplemental text B). We then obtained the differentialequation for ˆ h k ( t ) for k ≥ ddt ˆ h k ( t ) = ˆ f k − (cid:18) ρ sin πkN (cid:19) ˆ h k ( t ) + σ ˆ η k ( t ) . (7)ˆ h k ( t ) was explicitly derived by using Ito integral as fol-lows (See the supplemental text C):ˆ h k ( t ) = ˆ f k a k (1 − e − a k t ) − σ (cid:90) t e a ( s − t ) d ˆ B s . (8)ˆ B s is Brownian motion in the complex plane, and a k =4 ρ sin πkN . Equation (8) indicates that the Fourier coef-ficient ˆ h k ( t ) follows the complex Gaussian distributionwith the mean ˆ f k a k (1 − e a k t ) and the variance σ a k (1 − e − a k t ). The expected value of the power-spectrum | ˆ h k ( t ) | is written as E [ | ˆ h k ( t ) | ] = | ˆ f k | a k (1 − e − a k t + e − a k t )+ σ a k (1 − e − a k t ) . (9)Next, we introduce the squared local roughness w ( l, t )as the average of the variance of h x in m consecutive cells.Using the inter-cell distance ∆ x and non-negative integer m , w ( l, t ) is expressed as w ( l, t ) = w ( m ∆ x, t ) = R (0) − m ( m − m − (cid:88) l =1 m − l ) R ( l ) . (10) R ( l ) is an auto-correlation function of h x . Since R ( l ) isobtained by inverse Fourier transform of the power spec-trum | ˆ h k ( t ) | (Winner-Khinchin’s theorem), we obtainedthe expected value of w ( l, t ) as in (10). Figure 5(a)shows that the analytically derived w ( l, t ) was close tothe average of the numerical calculations of the circularmodel.Considering the slope of the log-log plots in Fig. 5(a),we obtained the Hurst and growth exponents as follows: α = log (cid:20) R (1) − R (2) R (0) − R (1) (cid:21) / β = t w ( l max , t ) ∂∂t w ( l max , t ) . (12)We calculated the expected values α = 0 .
80 and β =0 .
74 for the parameters used for the circular model inFig. 4. This confirmed that our analysis based on (6)captured the dynamics of the circular model in (3).Equation (11) also shows that the Hurst exponent isrepresented as the ratio of the linear sum of the powerspectra. Therefore, the multiplication of the power spec-trum by a constant value should not affect the Hurstexponents. By dividing the power spectra (9) by σ /ρ ,we found that the Hurst exponents were determined bythe index c and ρt , where c = ( v l − v f ) p l (1 − p l ) ρσ . (13)The plots of the Hurst exponents against log c with dif-ferent ρt are shown in Fig. 5(b). When the difference ofcell motility ( v l − v f ) is relatively large, c takes a large (b) (c)(a) FIG. 5. Parameter dependency of the Hurst and growth ex-ponents. (a) Comparison of numerically and analytically de-rived w ( l, t ). Log-log plot of w ( l, t ) against l at t = 18 (leftpanel), and of w ( l max , t ) against t (right panel). The light-colored lines indicate the numerically obtained values withdifferent sample paths, dots indicate the average of w ( l, t )for different paths, and the black line indicates analyticallyobtained values. (b) Numerically obtained parameter depen-dency of the Hurst exponent α on c and ρt . (c) The parameterdependency of the growth exponent β on ρct and ρt . Dots in-dicate the numerically obtained growth exponents, and theline indicates the plot of (14). value and α approaches α = 0 .
91, corresponding to thefixed noise model. On the other hand, when the randomcomponent of cell movement ( σ ) is relatively large, c issmall and α approaches α = 0 .
48, corresponding to EWmodel.The growth exponents β were also regarded as a two-variable function of ρt and c . From the theoretical consid-eration (See the supplemental text D), we can calculate β as follows: β = 2 (cid:0) − √ (cid:1) γ + 1 (cid:0) − √ (cid:1) γ + 4 . (14)where γ = ρct and we assumed that N and ρt are suffi-ciently large. Plots of β against log γ are shown in Fig.5( c ). It was shown that the random cell movement de-creases β . However, after sufficient time has passed, β always takes a value of 0 .
75, corresponding to fixed noisemodel.
E. Confirmation of analytical prediction
As described in this section, we experimentally exam-ined the dependency of the Hurst and growth exponentsas derived from numerical and analytical considerations.First, the effect of the initial size of the monolayer wasexamined. A monolayer of 3 mm diameter was prepared,and its spreading was observed. The Hurst and growthexponents were α = 0 .
837 and β = 0 . t = 4 to t = 18. The results of these observations are shownin Fig. 6(a) and Supplemental Movie 4. We foundthat the Hurst and growth exponents both decreased: α = 0 . , β = 0 . (cid:104) D ( θ ) (cid:105) was increased to 28 . µm , which is2 . α = 0 . ± .
012 and β = 0 . ± . γ = ρct such that the growth exponent β = 0 .
687 was satisfied in (14), and found γ = 17 . ρ and p l are regarded as identical to thosein the control (Fig. 4): ρ = 6 . , p l = 0 .
1. We obtained c = 0 .
148 and could then determine the values of σ and ( v l − v f ) p l (1 − p l ) to match the values of w ( l max , t ).In addition, the time evolution of the mean diameter is23 . µm/h , which is equal to v l p l + v f (1 − p l ), so both v l and v f can be estimated. The values for the estimatedparameters were v l = 123 , v f = 18 . σ = 31 . c = 0 .
19 from ρt = 120and α = 0 . v l = 126 , v f = 17 . σ = 28 . v l = 126 , v f = 17 . σ = 28 .
8. These estimated parameters reproducedthe experimental results.Next, we confirmed the consistency between theoreti-cally estimated and experimentally observed cell motil-ity. The experimental results (shown in Fig. 6(a)) werereproduced in our model by assuming that the motil-ity of the follower cell ( v f ) and the intensity of randomcell movement ( σ ) both increased about 9-fold from theparameters used in Fig. 2. Figure 6(c) and Supplemen-tal Movie 6 show the cell motilities on the different cellculture substrate. The total path lengths of the singlecell for 2 hr were measured: 6 . µm/h for the uncoated glass dish and 30 . µm/h for the collagen coated dish.The motility of the cells on the collagen increase to 4 . IV. DISCUSSION
In this study, we propose a concept that explains themechanism of the dynamic scaling law observed in thespreading of the epithelial cell monolayer. We found thatthe time evolution of the contour satisfied the dynamicscaling law: α = 0 .
86 and β = 0 .
73. Based on the ob-servation results, we constructed a simple mathematicalmodel, and demonstrated that the contour shape arosefrom the behavior of the cells at the sheet edge. Ourmathematical analysis of the model suggested that thepresence of leader cells is essential for observed patternformation. We found that the Hurst and growth expo-nents were dependent only on the ratio of the variance ofcell motilities to the intensity of the random cell move-ment. The theoretical prediction was experimentally con-firmed by changing the cell motilities. Thus, our studyoffers a new framework for examining dynamic scaling inthe biological phenomenon.The EW and KPZ models are widely known to satisfythe dynamic scaling law, and it has been reported thatthe contour shape of HeLa and Vero cell colonies followsKPZ universality. However, the MDCK cells observedhere did not follow this universality. Possible reasons forthis difference are the emergence of leader cells and thedifferences in observation time. First, the leader cellsemerge in the early stage and persist during observation,and then the differential motility corresponds to the fixednoise in the model. On the other hand, the EW and KPZequations do not contain the corresponding noise term.In addition, the effect of the nonlinear component in theKPZ system requires a long observation time. In thereport by Huergo and colleagues, the observation timewas 13000 minutes [20], while the observation time inthis study was 1080 minutes.We focused on the contour shape change that emergedfrom cell motility and cell-cell adhesion. In the view ofbiology, the model used in this study suggests that cellmotility played an important role in the contour forma-tion, while the effect of cell proliferation would be domi-nant after a long period of time. From a model perspec-tive, the increase in the number of cells on the circum-ference is not taken into account. Therefore, the modelcannot explain the pattern formation in a long period oftime that the effect of proliferation on the contour shape (c)(a)(b) V e l o c i t y ( μ m / h ) Control CollagenSlope = 0.687Slope = 0.742Slope = 28.7 (μm/hr) Slope = 0.689Slope = 0.741Slope = 28.3 (μm/hr) *** C on t r o l C o ll agen FIG. 6. The experimental result with 100 mg/L type I collagen coated dish is also reproduced by the model (2). (a)Experimental results. Time evolution of averaged D ( θ ) (left panel), log-log plot of the local roughness w ( l, t ) and l (middlepanel), and log-log plot of w ( l max , t ) and t (right panel). Dots indicate the experimentally obtained data and the dashed line isthe fitted line for the data. (b) Numerical results corresponding to (a). Parameters: v l = 126 [ µm h − ], v f = 17 .
87 [ µm h − ], p l = 0 . ρ = 6 . h − ], and σ = 28 . µm h − / ]. (c) Single cell movement assay. The trajectories of the cell nuclei areindicated by the colored lines in the uncoated glass dish (left panel) and the collagen coated dish (middle panel). Right panelshows the velocity of the cells; 6 . µm/h in the uncoated glass dish and 30 . µm/h in the collagen coated dish. *** p-value =3 . × − , scale bar = 100 µm . cannot be negligible.The growth exponents obtained from analysis of themathematical model, β = 0 .
74, were almost identical tothose obtained numerically and experimentally. Whilethe Hurst exponent, α = 0 . α = 0 .
85. The value of w ( l, t ) is defined as the mean of the standard deviationin the closed range l , however, this value cannot be di-rectly calculated from the power spectra. Therefore, toestimate the Hurst exponent, the expected value of vari-ance was calculated, and the square root of this valuewas taken to estimate the value of w ( l, t ). The valueof (cid:112) w ( l, t ) is not the same as w ( l, t ), since the meanof the standard deviation in a given interval is different from the square root of the mean of the variance. On theother hand, for the growth exponents, we calculate globalroughness w ( l max , t ) as the standard deviation with re-spect to the whole direction. Since this value is equal tothe square root of the variance, it is close to the valuesobtained numerically.The study of mathematical models of collective cellmovement has attracted significant attention over thepast decade, especially within the fields of statistical me-chanics and biophysics. The models can be categorizedinto continuous models and discrete models. In continu-ous models, the cell colony is regarded as a continuum.Such models have been constructed mainly to explain thefingering instability of the epithelial sheet [9, 28–32].0On the other hand, in discrete models, each cell hasbeen represented by polygons or particles. In the formercase, the classical vertex model with chemotaxis and fluidproperties [33] and the active vertex model [34], whichalso added the effect as an active fluid, have been pro-posed. In the latter, a particle model mostly included thecell-cell interactions and the random kinetic components.The model explicitly introduced leader cells [35], and amodel with the effects of the bending and the surfacetension have been proposed [36]. These discrete modelstend to be descriptive, and assume many factors that af-fect cell behavior. While such models are preferable toexplain the experimental data, however, it is likely to bedifficult to analytically explain the numerical results dueto the model complexity. Thus, analytical considerationsof the model equation are required to fully understandthe relationship between the physical quantity and thescaling property.In a sense, our model can be understood as the sim-plest form of the discrete particle model. The reason weused the spatially discrete model was to describe the dif-ferences of the motilities among the cells, such as leaderor follower cells. The result that even the simple modelexplained the dynamic scaling within the contour shapesuggested that the random cell movement, the determin-istic differences in cell motility, and the effects of intercel-lular adhesion have critical implications on the contourof the epithelial monolayer.The expected values of the power spectra (9) convergeto the constant values when ρt → ∞ . Since w ( l, t )are represented by the linear sum of the power spectra, w ( l, t ) also converges to constant values. Therefore, itis suggested that α at ρt → ∞ is determined by c , andthat there is a critical value c ∗ such that when c (cid:28) c ∗ ,the scaling law α = 0 . , β = 0 .
25 was obtained and when c (cid:29) c ∗ , α = 0 . , β = 0 .
75 was obtained. The value of c ∗ was estimated from the analytically obtained relationshipbetween expected values of w ( l, t ∞ ).The equation known as quenched Edwards-Wilkinson(QEW) equation [14, 37] is a model with the noise termdependent on x and h and the driving force F as follows: ∂∂t h ( x, t ) = ρ ∇ h ( x, t ) + F + η ( x, h ( x, t )) . (15)In this model, when F is sufficiently large, the noise termis considered to be spatially and temporally independent,then, the system could be identified as the EW equation: α = 0 . , β = 0 .
25. However, the transition to a differentdynamic scaling law α = 1 . β = 0 .
75 occur when
F (cid:28) F c [14, 37, 38]. The relationship between F c and thescaling law is similar to the relationship between c ∗ andthe scaling law in our model (6), which suggests that ourmodels could possibly be identified as QEW-type models.However, we expect that more complex theoretical meth-ods will be needed to solve this problem, and it remainsa topic for future research. ACKNOWLEDGMENTS
We acknowledge fruitful conversations with S. Ishihara(University of Tokyo), T. Ogawa (Meiji University) andH. Yasaki (Meiji University). This work has been sup-ported by JSPS through Grants No. JP18K06260 toH.T.I. [1] R. A. Weinberg,
The Biology of Cancer , 2nd ed. (GarlandScience, New York, 2013).[2] F. Lennon, G. Cianci, N. Cipriani, T. Hensing, H. Zhang,C. Chen, S. Murgu, E. Vokes, M. Vannier, and R. Salgia,Nat. Rev. Clin. Oncol. , 664 (2015).[3] O. Chepizhko, C. Giampietro, E. Mastrapasqua,M. Nourazar, M. Ascagni, M. Sugni, U. Fascio, L. Leggio,C. Malinverno, G. Scita, S. Santucci, M. Alava, S. Zap-peri, and C. Porta, Proc. National. Acad. Sci. , 11408(2016).[4] R. Mayor and E. Sandrine, Nat. Rev. Mol. Cell. Bio. ,97 (2016).[5] A. J. Ridley, M. A. Schwartz, K. Burridge, R. A. Firtel,M. H. Ginsberg, G. Borisy, T. J. Parsons, and A. Hor-witz, Science , 1704 (2003).[6] T. Omelchenko, J. Vasiliev, I. Gelfand, H. Feder, andE. Bonder, Proc. National. Acad. Sci. , 10788 (2003).[7] N. Yamaguchi, T. Mizutani, K. Kawabata, and H. Haga,Sci. Rep. , 7656 (2015).[8] J. Konen, E. Summerbell, B. Dwivedi, K. Galior, Y. Hou,L. Rusnak, A. Chen, J. Saltz, W. Zhou, L. Boise,P. Vertino, L. Cooper, K. Salaita, J. Kowalski, andA. Marcus, Nat. Commun. , 15078 (2017).[9] S. Mark, R. Shlomovitz, N. Gov, M. Poujade, G. Erwan, and P. Silberzan, Biophys. J. , 361 (2010).[10] M. Vishwakarma, J. Russo, D. Probst, U. Schwarz,T. Das, and J. Spatz, Nat. Commun. , 3469 (2018).[11] R. Riahi, J. Sun, S. Wang, M. Long, D. Zhang, andP. Wong, Nat. Commun. , 6556 (2015).[12] R. F. Voss, Physica D , 362 (1989).[13] B. B. Mandelbrot, Phys. Scr. , 257 (1985).[14] P. Meakin, Fractals, Scaling and Growth Far from Equi-librium (Cambridge University Press, Cambridge, 2011).[15] J. Moreira, J. da Silva, and S. Kamphorst, J. Phys.Math. Gen. , 8079 (1994).[16] J.-i. Wakita, H. Itoh, T. Matsuyama, and M. Matsushita,J. Phys. Soc. Jpn. , 67 (1997).[17] S. Santalla, R. Javier, J. Abad, I. Mar´ın, M. Espinosa,M. Javier, L. V´azquez, and R. Cuerno, Phys. Rev. E ,012407 (2018).[18] M. Myllys, J. Maunuksela, M. Alava, A. T, J. Merikoski,and J. Timonen, Phys. Rev. E , 036101 (2001).[19] K. Takeuchi and M. Sano, J. Stat. Phys. , 853 (2012).[20] M. Huergo, M. Pasquale, P. Gonz´alez, A. Bolz´an, andA. Arvia, Phys. Rev. E , 011918 (2012).[21] N. Muzzio, M. Pasquale, P. Gonz´alez, and A. Arvia, J.Biol. Phys. , 285 (2014).[22] S. F. Edwards and D. Wilkinson, Proc. Royal. Soc. Lond. Math. Phys. Sci. , 17 (1982).[23] M. Kardar, G. Parisi, and Y. Zhang, Phys. Rev. Lett. , 889 (1986).[24] D. Roy and R. Pandit, Phys. Rev. E , 030103 (2019).[25] K. A. Takeuchi, Physica A , 77 (2018).[26] J. Tinevez, N. Perry, J. Schindelin, G. Hoopes,G. Reynolds, E. Laplantine, S. Bednarek, S. Shorte, andK. Eliceiri, Methods , 80 (2017).[27] H. Haga, C. Irahara, R. Kobayashi, T. Nakagaki, andK. Kawabata, Biophys. J. , 2250 (2005).[28] P. Carlos, R. Alert, B. Carles, G. Manuel, T. Kolodziej,E. Bazellieres, J. Casademunt, and X. Trepat, Nat. Phys. , 79 (2019).[29] G. Ouaknin and B. Pinhas, Biophys. J. , 1811 (2009).[30] P. Lee and C. Wolgemuth, Plos. Comput. Biol. ,e1002007 (2011). [31] M. H. K¨opf and L. M. Pismen, Soft Matter , 3727(2013).[32] R. Alert, B. Carles, and J. Casademunt, Phys. Rev. Lett. , 088104 (2019).[33] M. Salm and L. Pismen, Phys. Biol. , 026009 (2012).[34] D. L. Barton, S. Henkes, C. J. Weijer, and R. Sknepnek,Plos. Comput. Biol. , e1005569 (2017).[35] N. Sep´ulveda, L. Petitjean, O. Cochet, G. Erwan, P. Sil-berzan, and V. Hakim, Plos. Comput. Biol. , e1002944(2013).[36] V. Tarle, A. Ravasio, V. Hakim, and N. S. Gov, Integr.Biol. , 1218 (2015).[37] D. A. Kessler, H. Levine, and Y. Tu, Phys. Rev. A ,4551 (1991).[38] O. Narayan and D. S. Fisher, Phys. Rev. B48