Choice of friction coefficient deeply affects tissue behaviour in stochastic epithelial vertex models

To understand the mechanisms that coordinate the formation of biological tissues, the use of numerical implementations is necessary. The complexity of such models involves many assumptions and parameter choices that result in unpredictable consequences, obstructing the comparison with experimental data. Here, we focus on vertex models, a family of spatial models used extensively to simulate the dynamics of epithelial tissues. Usually, in the literature, the choice of the friction coefficient is not addressed using quasi-static deformation arguments that generally do not apply to realistic scenarios. In this manuscript, we discuss the role that the choice of friction coefficient has on the relaxation times and consequently in the conditions of cell cycle progression and division. We explore the effects that these changes have on the morphology, growth rate and topological transitions of the tissue dynamics. These results provide a deeper understanding of the role that an accurate mechanical description plays in the use of vertex models as inference tools. This article is part of a discussion meeting issue ‘Causes and consequences of stochastic processes in development and disease’.


Introduction
In order to model the dynamics of epithelial tissues, the use of vertex models is increasing.Successful inference from such models is challenging due to the large number of implicit parameters.This has resulted in a focus of the scientific community on studying a subset of the parameters of vertex models while ignoring the effects of others.In this article, we focus on exploring the role of one of these neglected parameters, the friction coefficient.This coefficient has been traditionally neglected using arguments of quasi-static deformation processes, which do not account for tissue viscosities.In the analysis performed in this paper, we show computationally and analytically how this coefficient's choice impacts the tissue's morphology, the cell cycle dynamics, and the topology transitions of the tissue.This provides a deeper understanding of the role that an accurate mechanic description plays in the use of vertex models as inference tools.
The spatio-temporal organization of epithelial cells is central to many biological processes.During embryo development, the formation and growth of epithelial layers orchestrate the axes formation and establish the cell lineages in the embryo [1][2][3][4][5][6][7][8][9][10][11][12][13].In addition, the mechanics of epithelial structures play a major role in tumour growth and malignancy [14].The geometrical simplicity of epithelia has allowed the development of many mathematical models to understand the tissue's morphology.Capturing the relationship between the biophysical properties of individual cells and tissue-level properties [10,[15][16][17][18], such an approach provides a powerful tool to acquire generalistic insight into the tissue.
Unfortunately, direct application of analytical frameworks to particular biological tissues is often not possible.The complexity resulting from the highly dynamic biological processes involved, such as the coupling of cell mechanics and cell cycle progression, makes necessary the development of numerical tools capable of recapitulating the growth of the tissue [7].In this aspect, computational models have become an indispensable tool to predict and distinguish between hypothesized biophysical mechanisms, as well as to plan future experiments.One particular family of approaches that have proven to be highly successful to model epithelial tissues is that of vertex models [5,[19][20][21][22].These models assume packing geometries in which cell sheets are approximated by tessellations of polygons [2,16,23,24].Under this prescription, cells can be fully described as a set of vertices.Tissue dynamics can be naturally incorporated by formulating the evolution of the vertices of the cells as an energy minimization problem [16] resulting in the relaxation of a prescribed Hamiltonian.The flexibility of this approach allows us to include tailored biophysical energetic terms accounting for biophysical properties of the cells such as adhesion and contractibility.Additionally, it also enables the inclusion of cell-specific physiological processes such as cell-cycle progression, differentiation or apoptosis.
Traditional approaches describe the tissue morphology through the quasi-steady state solution of the Hamiltonian [15,16,25].Nevertheless, this quasi-steady solution does not allow the inclusion of active deformations of the tissue.Thus, a formulation of the equations of motion dictating the evolution of the tissue is necessary.The most adopted solution to this problem is to assume that all the vertices follow overdamped dynamics.Collectively, cells passively respond to the resulting forces in a viscoelastic manner.Hence, to reveal the mechanics of tissue morphogenesis, we must understand the relationship between active force generation and passive viscoelastic response at the multicellular level.While this bestows the model with the ability to reproduce more realistic dynamics, the effect that different choices of viscosity have on the resulting dynamics is rarely explored [26].This choice of friction will not be relevant in tissues where the mechanic relaxation timescales are faster than the rest of cellular processes (as shown in ablation experiments in the Drosophila wing disc [16]).On the other hand, choice of friction might not be overlooked in situations where mechanical timescales are comparable with those of the cell-cycle dynamics resulting in a continuous nonlinear coordinated effect of both processes [7,12,13].In this manuscript, we explore the effects of the choice of friction, emphasizing its effects in tissue evolution such as variability of cell shape, stochastic cell-cycle duration, tissue growth and topological cell rearrangements.

(a) Vertex dynamics
In order to describe the evolutionary dynamics of the characteristic polygonal morphology of a planar vertex model, the forces induced by the prescribed energy potential can be encoded into an equation of motion for the N vertices describing the tissue r ¼ fr i jr i [ R 2 , i ¼ 1, . .., Ng.The dynamics of every single vertex can be described as a deterministic overdamped equation of motion [25], where inertial terms are neglected compared to dissipative terms, and the friction parameter γ is constant for every single vertex in the system, where F i denotes the total force acting on vertex i at time t, and it can be directly derived from the energy potential E as F i ¼ Àr i E (figure 1a).The choice for the energy potential includes the biophysical properties of the cell relevant for the particular tissue under study and usually has the form [16], where the first sum runs over all the cells α in the tissue incorporating cell-intrinsic energy terms of elasticity and contractility.The elasticity depends on the deviation from an autonomous preferred area A 0 a ðtÞ modulated by an elasticity constant K a .Similarly, the contractility of the cell is introduced as an elastic term that depends on the perimeter of the cell L a and a contractility constant G a .The second sum in equation (1.2) runs over all the edges 〈ij〉 incorporating cell interaction terms such as cell-cell adhesion energy, where L ij is a constant representing the line tension.More sophisticated expressions for the potential can include other effects such as inhomogeneities of the contractile tension or non-harmonic energy terms.

(a) Effect of friction coefficient on cell morphology
The prescribed equations of motion given by equations (1.1) and (1.2) allow a reparametrization of time τ = t/γ that suggests that the role of the friction coefficient γ does not perturb the expected configuration of the tissue.Under this royalsocietypublishing.org/journal/rstbPhil.Trans.R. Soc.B 379: 20230051 assumption, the analysis of the morphology of the tissue is usually reduced to understanding the roles of adhesion and contractility [2,16,24,27].Unfortunately, this dynamical symmetry breaks when other timescales influence the behaviour of the tissue.The main mechanism by which this happens is cell proliferation in which cells divide and change their morphology as they progress through the cell cycle.This is included in vertex models through a time dependence in the energy term (usually encoded by a change in the optimal cell area A 0 (t)).This active deformation of the cell not only will directly affect its shape but will also affect the topology of the network when cell divisions occur.How relevant these effects are in the epithelium configuration will depend on the relationship between the different timescales.
In addition to cell-autonomous timescales, it is common to assume that there are physical constraints to cell-cycle progression such as the requirement for the cell to reach a certain critical size, A c , and a minimum cell-cycle time of average t c for cytokinesis to take place [12,16,[27][28][29][30].All these factors introduce feedback between the timescales of relaxation and cell-cycle progression obscuring the effect that the choice of γ has on the epithelium morphology.In order to explore these different effects we made use of the vertex model implementation TifoSi [30].Specifically, in order to unravel the mechanic's feedback on the cell-cycle progression we used the simplified prescribed target area, where t c is the intrinsic cell-cycle length of the cell.To incorporate variability in cell-cycle progression among different cells, t c is described as a stochastic variable assigned to each cell at birth.We described t c as a sum of a deterministic minimum cell-cycle duration and a stochastic exponential event [30], with T c the average threshold time required for a cell to divide (〈t c 〉 = T c ), and 1 [ ½0, 1 a parameter that weights the contribution of the deterministic and the stochastic part.Since t c is assigned to each cell independently of the interaction with other cells, we will refer to this variability as intrinsic.A summary of the values used throughout the manuscript can be found in table 1, other details can be found in the electronic supplementary material.
To examine the potential effects that the choice of friction has on the tissue morphology, we compared the resulting epithelia for different values of the friction constant γ that start with an identical hexagonal lattice of 100 cells for a duration of t = 3200 ≃ 10 T c .Visual inspection of the resulting tissues reveals that cell size variation increases with the friction coefficient (see figure 1b).
Further quantitative analysis of these differences shows that despite area heterogeneities, the tissues are geometrically very similar.The distribution of relative areas as a function of a number of neighbours follows the quadratic property of the cell arrangement as has been described computationally and experimentally in Kokic et al. [31].Discrepancies only appear for larger cells with more than seven neighbours (see figure 2a).Similarly, the average shape index of the cells defined as the ratio between the apical perimeter of the cell and the square root of its area remained independent of the friction value chosen, suggesting that the measured change in cell area across friction values is isotropic (see figure 2b).In addition, the shape index has been identified in previous studies with the change of phase in the tissue [32,33].For the parameters explored, we obtained values below the solid-liquid transition, suggesting that the friction of the tissue did not affect the phase behaviour of the tissue that stayed at every moment in a solid (jammed) phase.Exploration of other parameters such as distribution of polygon shapes, cell area at fixed final tissue size, and cell perimeter distributions depict a similar scenario (electronic supplementary material, figures S1-S3).All in all, the morphology of the cells seems to be conserved for different values of friction, with the main difference being an overall slight decrease in the size of the cells for increasing values of friction.Specifically, we observed a reduction of 20% of the average apical area over the three orders of magnitudes of values of friction explored (figure 3a).

(b) Effect of friction coefficient on tissue growth
The reduction in apical cell area as friction increases can be understood through the relationship between friction and apical cell junction rearrangement required for mitosis.This was confirmed upon inspection of the apical areas of the cells at the moment of division (figure 3c).Furthermore, analysis of the age of the cell upon division reveals that cells with lower friction values divide faster than high friction ones (figure 3b).This inverse relationship between size and age of division emerges from the prescribed conditions for cytokinesis.Lower frictions allow for a fast rearrangement of the tissue allowing mitotic cells to reach the minimum required division area A c before the cell-cycle minimum duration t c .On the other hand, for higher frictions, the rearrangement is slower and cells require time beyond t c to reach the apical cell division area A c .Thus, the tissue friction γ controls the interplay between spatial and temporal cytokinesis conditions, attaining the limiting temporal or spatial scenarios for, respectively, low or high values of γ.
This effect that γ has on the cell cycle results in a change in the average cell-cycle length and consequently a dramatic change in the tissue size at different times (figure 4).In Table 1.Simulation parameter values.If not stated explicitly, the values are the same as the ones in Canela-Xandri et al. [30].where N 0 is the initial number of cells and 〈T〉 cell is the average cell-cycle duration obtained in the simulations.The third expression in equation (2.3) is identical to the second but expressing the growth in terms of the growth rate k det ≡ ln(2)/ 〈T〉 cell .Strikingly, the deterministic growth (equation (2.3)) predicts a growth faster than the computational observation (dashed lines in figure 4).We hypothesized that not only the change in the mean cell-cycle length but also the spread of the distribution of cell-cycle durations have a relevant effect on the cell-cycle duration and cannot be neglected.For a stationary cell-cycle distribution the growth of the tissue will still be exponential but with an effective rate k sto , NðtÞ ¼ N 0 e kstot : ð2:4Þ In order to relate k sto with the distribution of cell-cycle durations we can write a continuity equation for the density of cells N at time t with an age τ, where the age of a cell is the time elapsed since the birth of that cell [34].The density N ðt, tÞ is related to the cell number N (t) through normalization, In a small window of time δ the density of cells N ðt, tÞ is reduced an amount proportional to the probability g d ðtÞ that a cell of age τ divides in that time interval, ð2:6Þ The probability g δ (τ) of division at an age τ, can be written in terms of the distribution of cell cycle durations P(τ) that does not change as the tissue grows, similar to the distributions observed in the vertex model simulations.In order to do so, we can write g δ (τ) as where GðtÞ ¼ Ð t 0 PðtÞ dt is the cumulative density function of cell cycle durations.Introducing the expression for the division probability (equation (2.7)) in the continuity equation (2.6) and taking the limit δ → 0, we get the differential form of the continuity equation, The first term of equation (2.8) can be related to the growth rate of the tissue k sto ¼ ð1=N ðt, tÞÞð@N ðt, tÞ=@tÞ.This gives rise to an ordinary differential equation for the age τ,  Note that there is not an immediate relationship between k sto and k det , the latter depending exclusively on the first moment of the cell-cycle distribution hTi cell ¼ Ð 1 0 tPðtÞ dt.For the intrinsic cell cycle distribution (equation (2.2)), in the absence of cellcell interaction (〈T〉 cell = T c ), equation (2.12) becomes, In the deterministic cell-cycle scenario (1 ¼ 1), we recover the deterministic growth k sto = ln2/T c .On the other hand, for an exponentially distributed cell-cycle duration (1 ¼ 0) we obtain an upper bound for the tissue growth rate k sto = 1/T c .For intermediate values of 1, the rate k sto can be obtained by solving numerically the transcendental equation (2.13).For the parameters used in this paper (table 1), we have k sto T c ≃ 0.70, which predicts a tissue growth very close to the deterministic scenario.In general, for the stochastic cell-cycle model used equation (2.2) it is expected that the intrinsic stochasticity will increase the value of k sto , i.e. while preserving the same average cell cycle duration, the tissue will grow at a faster speed than in a purely deterministic scenario, which is the opposite effect to the slow down predicted by the simulation of the vertex model (solid line in figure 4).This apparent mismatch is not surprising, the distribution of cell-cycle durations in the simulation will be different from the intrinsic stochastic cell cycle in the absence of cell-cell interactions (cf.equation (2.2) and figure 3b).Since we cannot access analytically the expression for cell-cycle durations in the tissue P(τ), the actual value k sto can be calculated by solving numerically the integral equation (2.12).Alternatively, one can integrate the continuity equation (2.8) for all the possible ages τ, and use the normalization of N and the proliferation boundary condition equation (2.11) obtaining the relationship, that relates directly the rate of growth of the tissue and the age structure of the population.Comparison of the rates of growth of the tissue with the computational observation is able to recover the discrepancies that appeared with the analysis (dashed line in figure 4).This reveals the relevance that cytokinesis conditions impose on the growth of the tissue by shaping the cell-cycle length distribution.
(c) Topological transitions are also affected by the friction constant In the previous section, we showed that friction has a dramatic effect on the growth of the tissue by affecting the prescribed typical mitosis conditions used in the vertex model.In particular, the threshold condition in critical cell area A c seems to have a relevant impact on the mean duration of the cell cycle and its variability.To explore further this effect we explored the resulting tissues from a model that does not require a threshold in the target area in order to divide (A c = 0).As expected, the resulting tissue grows faster, following a growth more similar to the intrinsic stochastic growth model equation (2.13) (figure 5a).Nevertheless, the agreement with the deterministic model is not perfect, suggesting that more factors affect the timely mitosis of the cells.Analysis of the configuration of cell areas confirms that removing the condition of critical area affects the distribution of cell sizes for large values of γ but does not rescue the distribution of cell sizes observed at lower values of γ (figure 5b).Cells do not require any more to reach a certain size to divide, which results in overall smaller cell sizes at division, and consequently smaller cell areas at the start of the cell cycle of both daughter cells.Since, in the absence of a cell critical area, cells progress along the cell cycle without affecting each other, the cell area should not affect the timing of division and tissue growth.The only mechanism left by which cells will change their cell-cycle progression is a topological transition.Specifically, T2 transitions control the removal of cells when they become too small.This mechanism has been employed in the past to explain apoptosis or extrusion of cells from the epithelium [12].Analysis of the amount of T2 transitions confirms that the removal of A c results in an increased number of T2 transitions sufficient to explain the discrepancies between simulations with different friction coefficients (figure 5).This is only observed for high values of γ where the slow vertex dynamics impede the growth of the cells up to the target area A 0 .

Discussion
Vertex models are a very powerful computational tool to simulate epithelial organization.One of their advantages is the possibility of incorporating intrinsic mechanistic details of cell-cycle progression through a target area term in the Hamiltonian [12,30].This allows us to understand the nonlinear effects by which cell-intrinsic mechanistic details affect overall tissue dynamics.At the same time, it reveals how sensitive the results are to different modelling choices that are usually taken arbitrarily in the literature.
To make this evident, in this paper we have focused on the effects of the choice of friction coefficient γ in the evolution of the tissue (equation (1.1)).In particular, we have observed that not only do different frictions result in nonintuitive changes in size and variability of cell areas, but most importantly in tissue growth.Typical prescriptions of the cell cycle in vertex models are limited to incorporating average cell observations in the form of a time-dependent target area and a set of conditions of division.By contrast here we show that the choice of friction will affect dramatically the resulting cell cycle by controlling these division conditions and modifying its mean duration and dispersion.All in all, this reveals the complexity of incorporating intrinsic cell-autonomous mechanistic descriptions into sophisticated vertex models and how accurate descriptions of the variability of cell-cycle duration in time are paramount to achieving a useful description of the tissue.This is in line with previous results [35], where they explored the time-scale separation by varying cell-cycle duration in a non-dimensional model while keeping the friction constant, observing a substantial impact in tissue configuration.In this particular case, they royalsocietypublishing.org/journal/rstbPhil.Trans.R. Soc.B 379: 20230051 also observed an impact on the fluidity of the tissue through the increase of topological transitions.These results are of particular relevance given recent observations on embryo development where fluidity transitions are thought to control cell organization in the mouse neuroepithelium [13] and in vertebrate body axis elongation [36].
Further efforts into a proper description of the cell cycle need to incorporate the stochastic biological mechanisms of cell-cycle progression.While in this study we followed custom stochastic descriptions for minimum threshold times of division, different models are compatible with the intrinsic stochasticity in cellcycle advance [34].Actual applications to specific tissues will need a better understanding of the distributions of cell-cycle duration at a molecular level and interrogate carefully how these are related to mechanical feedbacks.
As part of the integration of these nonlinear effects in tissue growth, we showed that stochastic distributions of the cell cycle can be tackled through continuity equations for the age of the cells in the tissue.This approach can also accommodate the effect of T2 transitions by incorporating cell clearance along its cell cycle.For instance, an asymmetric division in which only a fraction of daughter cells stays in the tissue after division [37] could be incorporated replacing the factor 2 in equation (2.11).All in all, these calculations highlight the relevance of measuring full distributions of cell-cycle lengths in experimental assays as opposed to reporting summary statistics.
The current study has tackled the usual dissipative prescription by which the displacement per unit of time is proportional to the force homogeneously in the tissue (equation (1.1)) where γ is a global constant of the tissue.Future studies should also question this approximation that is based on the implicit physical assumption that each vertex is a point of mass with identical assigned friction.By contrast, one can expect that the dissipation of each vertex is related to the total material that needs to be displaced under a certain force, e.g. it is not the same to drag the membranes joining three small cells, than the membranes of three bigger ones.
All in all, the ultimate goal of vertex models is to be used together with experimental data to understand the mechanistic properties of tissue dynamics.The findings of this paper, combined with recent research highlighting a substantial zone of non-identifiable parameters in vertex models compatible with a large variety of biological scenarios [38], underscore the urgent need for further work in this field.This effort is essential to reliably use vertex models as tools for inferring mechanistic properties from experimental data.royalsocietypublishing.org/journal/rstb Phil.Trans.R. Soc.B 379: 20230051

Material and methods
All the simulations were produced using the vertex model implementation TiFoSi.An example of the TiFoSi input .xmlfile used to run the simulations is attached as electronic supplementary material.The input .xmlfile contains all the parameters required to run the simulations with the most important parameters summarized in table 1 and equations (2.1), (2.2).

Figure 1 .
Figure 1.(a) Force scheme applied in each vertex of a cell.(b) Tissue plot for different values of the friction γ : 0.01, 0.1, 1 and 2. Simulation duration t = 1600 ≃ 5 T c .Other parameter values are provided in table 1.

Figure 2 .
Figure 2. (a) Average cell area A n normalized by the mean area A as a function of the number of neighbours, n.(b) Average shape index, s n ¼ P n = ffiffiffiffi A n p , with P the perimeter of the cell, as a function the number of neighbours, n.Each line is the average of 20 different realizations with different values of γ : 0.01 (red), 0.1 (purple), 1 (blue) and 2 (green).Error bars indicate the standard error of the mean (s.e.m.).Simulation duration t = 3200.Other parameter values are provided in table 1.

Figure 3 .
Figure 3. Variation of cell size and division properties.(a) Distribution of cell areas in the tissue aggregated for 20 different simulations at t = 3200.(b) Division time distribution for different values of γ at different tissue sizes.Lines indicate the median (solid line) and quartiles (dotted lines) for the distribution resulting from aggregating 20 different realizations.(c) Comparison of cell area and age at division time.1T c is the minimum duration of any cell cycle (see equation (2.2)) and A c is the critical area required to divide.Simulation corresponds to 20 tissues of 15 000 cells aggregated.Parameter values are provided in table 1.

Figure 4 .
Figure 4. Comparison of the simulated tissue size (circles) with different deterministic and stochastic predictions (lines) for different values of the friction coefficient.The stochastic intrinsic model (solid line) corresponds to the model without cell-cell interaction, equation (2.13).The deterministic tissue model (dotted line) corresponds with an exponential growth model with duplication time obtained from the average in the tissue (〈T〉 cell ), equation (2.3).The stochastic tissue prediction (dashed line) was obtained by introducing the numerical age structure at the last snapshot of the tissue in equation(2.14).Tissues were simulated until a size of 15 000 cells.Other parameter values are provided in table 1. Simulated circles are the average of 20 different realizations, standard error of the mean is smaller than the circle size.

Figure 5 .
Figure 5.Effect of topological transitions in tissue growth.(a) Comparison of tissue growth with and without critical area conditions for mitosis (dotted and solid lines) for γ = 1.Both conditions result in a tissue smaller than expected from the intrinsic growth (dashed line) resulting from equation (2.13).Inset: zoom on the diverging parts of the curves using a linear scale.Results are the average of 20 simulations for each condition.(b) Distribution of cell areas in the tissue at 15 000 cells for different values of γ and condition of critical division area.Violin plots show the aggregated data from 20 simulations for each condition.(c) Number of T2 transition events in the absence of critical area condition for two different values of γ.Each condition is the aggregate of 20 realizations (circles).All violin plots indicate the median (solid line) and the quartiles (dotted lines).Other parameter values are provided in table 1.
).The deterministic tissue model (dotted line) corresponds with an exponential growth model with duplication time obtained from the average in the tissue (〈T〉 cell ), equation(2.3).The stochastic tissue prediction (dashed line) was obtained by introducing the numerical age structure at the last snapshot of the tissue in equation (2.14).Tissues were simulated until a size of 15 000 cells.Other parameter values are provided in table 1. Simulated circles are the average of 20 different realizations, standard error of the mean is smaller than the circle size.royalsocietypublishing.org/journal/rstb Phil.Trans.R. Soc.B 379: 20230051