Impact of implementation choices on quantitative predictions of cell-based computational models

‘Cell-based’ models provide a powerful computational tool for studying the mechanisms underlying the growth and dynamics of biological tissues in health and disease. An increasing amount of quantitative data with cellular resolution has paved the way for the quantitative parameterisation and validation of such models. However, the numerical implementation of cell-based models remains challenging, and little work has been done to understand to what extent implementation choices may influence model predictions. Here, we consider the numerical implementation of a popular class of cell-based models called vertex models, which are often used to study epithelial tissues. In two-dimensional vertex models, a tissue is approximated as a tessellation of polygons and the vertices of these polygons move due to mechanical forces originating from the cells. Such models have been used extensively to study the mechanical regulation of tissue topology in the literature. Here, we analyse how the model predictions may be affected by numerical parameters, such as the size of the time step, and non-physical model parameters, such as length thresholds for cell rearrangement. We find that vertex positions and summary statistics are sensitive to several of these implementation parameters. For example, the predicted tissue size decreases with decreasing cell cycle durations, and cell rearrangement may be suppressed by large time steps. These findings are counter-intuitive and illustrate that model predictions need to be thoroughly analysed and implementation details carefully considered when applying cell-based computational models in a quantitative setting.

developmental biology, the larval wing disc of the fruit fly Drosophila [3, 4, 25]. We conduct 87 convergence analyses of vertex positions with respect to all numerical and non-physical model 88 parameters, and further analyse to what extent experimentally measurable summary statistics 89 of tissue morphology, such as distributions of cell neighbour numbers and areas, depend on 90 these parameters. 91 We find that vertex model predictions are sensitive to the length of cell cycle duration, the 92 time step, and the size of the edge length threshold for cell rearrangement. Specifically, vertex 93 configurations do not converge as the time step, the edge length threshold for cell rearrangement, influence the rate of cell rearrangement. Counterintuitively, the rate of cell removal is robust 98 to changes in the area threshold for cell removal over multiple orders of magnitude. Further, 99 analysing the active forces within the tissue reveals that vertices are subject to stronger forces 100 during periods when cells grow and divide. 101 The remainder of the paper is organised as follows. In section 2, we describe our vertex where µ is the friction strength, x i (t) is the position vector of vertex i at time t, and E denotes 117 the total stored energy. The number of vertices in the system may change over time due to 118 cell division and removal. The symbol ∇ i denotes the gradient operator with respect to the 119 coordinates of vertex i. The total stored energy takes the form is its target area. This term penalises deviations from the target area for individual cells, thus 122 describing cellular bulk elasticity. The second sum runs over all cell edges i, j in the sheet and 123 penalizes long edges (we choose Λ > 0), representing the combined effect of binding energy and 124 contractile molecules at the interface between two cells. The third sum also runs over all cells, 125 and P α denotes the perimeter of cell α. This term represents a contractile acto-myosin cable (1) and (2) become where x i , A α , A 0,α , l i,j and P α denote the rescaled i th vertex positions, the rescaled area and  To solve equations (3) and (4) numerically we use a forward Euler scheme: We analyse the dependence of simulation outcomes on the size of ∆t in the Results section.  The assigning of these cell cycle stages to two thirds and one third of the total cell cycle 181 duration t l , respectively, allows us to modify the average age of a dividing cell with a single

211
In this section, we analyse how model behaviour depends on numerical and non-physical model

228
We analyse the dependence of the summary statistics on the cell cycle duration, t l , in figure 2.

229
The cell number and tissue area at the end of the simulation, and the total number of cell that each cell has. This summary statistic is often used to characterise epithelia [3, 26, 54,55]. 259 We find that the mean cell area for each polygon number is not sensitive to changes in cell cycle 260 length and increases monotonically with polygon number (figure 2F). 261 We interpret the data in figure 2 as follows. Differences in tissue size and cell packing arise    x (t + ∆t ) = x (t ) + ∆t 6 (k 1 + 2k 2 + 2k 3 + k 4 ) , Here, ∇ denotes the gradient with respect to the vector x.     prohibitive increases in calculation times as the time step is decreased. In previous studies, 506 vertex models have been reported to converge as the time step is decreased [45,56]. Our 507 analysis differs from these previous studies by considering a tissue undergoing cell division and 508 rearrangement rather than relaxation from an initial condition.

509
The simulation results are sensitive to the T1 transition threshold chosen in the simulation.

510
The size of the T1 transition threshold can be used to regulate the extent to which the simulated 511 tissue is allowed to rearrange in order to minimise energy. Literature values for this quantity