Quantitative cell-based model predicts mechanical stress response of growing tumor spheroids over various growth conditions and cell lines

Model simulations indicate that the response of growing cell populations on mechanical stress follows the same functional relationship and is predictable over different cell lines and growth conditions despite experimental response curves look largely different. We develop a hybrid model strategy in which cells are represented by coarse-grained individual units calibrated with a high resolution cell model and parameterized by measurable biophysical and cell-biological parameters. Cell cycle progression in our model is controlled by volumetric strain, the latter being derived from a bio-mechanical relation between applied pressure and cell compressibility. After parameter calibration from experiments with mouse colon carcinoma cells growing against the resistance of an elastic alginate capsule, the model adequately predicts the growth curve in i) soft and rigid capsules, ii) in different experimental conditions where the mechanical stress is generated by osmosis via a high molecular weight dextran solution, and iii) for other cell types with different growth kinetics from the growth kinetics in absence of external stress. Our model simulation results suggest a generic, even quantitatively same, growth response of cell populations upon externally applied mechanical stress, as it can be quantitatively predicted using the same growth progression function.


Center-based model: friction terms and forces
We look at the Equation of motion (4) for cells in more detail. The general form for the friction tensors in Equation (4) reads with u ij = ( r j − r i )/|| r j − r i ||, where r i , r j denote the position of the centers of cell i and object j. I is the 3 × 3 identity matrix, ⊗ denotes the dyadic product [12]. The individual cell friction coefficients are γ ⊥ and γ || , respectively perpendicular and parallel to the movement direction. As experiments do not indicate ECM inhomogeneity or anisotropy, the cell-ECM friction matrix is considered to be diagonal. The ECM was not represented explicitly as experimental observations indicated that the ECM fraction is small and approximately homogeneously and isotropically distributed in the intercellular spaces. As represented in Equation 5 the interaction force resulting from compression, deformation and adhesion can be expressed as a sum of a repulsion force, here represented by a modified Hertz contact force, and an adhesion force. The interaction force acts along the line connecting the centers of two interacting spherical cells.
The micro-motility force of cell i is mimicked by a Brownian motion term with zero mean value and are uncorrelated in time: < F mig,i >= 0 with < F (t 1 ) mig,i ⊗ F (t 2 ) mig,i >= Aδ(t 1 − t 2 ). (0.2) Based on a formal analogy to colloidal particle systems [19], the autocorrelation amplitude of the force as approximated A = 2Dγ 2 I, D being the cell diffusion constant, γ the friction coefficient of a cell in the medium and I the unity matrix. The scalar γ emerges in an isotropic environment, for which γ || = γ ⊥ (:= γ). In an extracellular matrix environment, A is largely controlled by the cell itself [3]. The equations of motion (4) do not conserve total momentum due to the micro-motility term as part of the momentum is transferred to the ECM, that is not explicitly modeled here. The system of Equation 4 was integrated numerically until the simulation time surpasses the total duration of the experiment which is ∼ 10d. Equation 4 results in a linear problem with a sparse symmetric positive definite matrix, which can be solved efficiently by a Conjugate Gradient method [12,13] with Jacobi preconditioner to obtain the cell velocities and an explicit Euler integration scheme to obtain the cell postions.

Deformable Cell Model: friction terms and forces
We recapitulate the basic equation of motion for a DCM denoted in (Equation 13). The individual elements generating the in-plane elastic forces between the surface nodes represented by the 1st terms on the lhs and rhs of Equation 13 are modeled by classical linear spring-damper systems. The force between the nodes captures the elastic response of the shell-like structure including the cortical cytoskeleton of the cell. When the elastic and dissipative components are summed up, one acquires a Kelvin-Voigt element. The vector force between two nodes i and j reads (see Figure A): where F e,ij and F v,ij are the elastic and dissipative forces, k s is the spring constant, γ represents the dissipation, d 0,ij , d ij are the initial and actual distance vectors, and v ij is the relative velocity between the nodes, respectively. The second term in Equation 13 representing the surface bending resistance is incorporated by the rotational resistance of the hinges determined by two adjacent triangles α = ijk and β = ijl . This defines the bending moment M : where k b is the bending constant, and θ is determined by the normals to the triangles n α , n β , with θ 0 being the angle of spontaneous curvature. M can be transformed to an equivalent force system F m acting on each of the nodes of the triangles [13]. Restoring volume compression / expansion forces F vol,i controlled by the bulk compression modulus of the cytoplasm are computed from the internal pressure using the volume change and bulk modulus of the cytoplasm K. In our simulations K reflects the overall cell compressibility including the cortex, as in its physiological range of elasticity, the cortex contributes little to the overall bulk modulus of the cell, i.e. K E cor h cor /R cell ( [8], own test runs). In analogy to Equation 8 the pressure applied to the cell is therefore approximately given by: The forces F vol on the nodes are perpendicular to the cell surface and can be obtained by multiplying the pressure with the surface area assigned to each node (each node has 6 adjoining triangles which each count for 1/3 of the surface area).
During the simulations, large variations in area of the triangles in the network can cause numerical artifacts. These can be avoided by adding a force F T , which is proportional to the area expansion of the individual triangles [15,20]: Here, A 0 and A are the initial and the current areas of the triangle, and k A is the area compression stiffness. The forces F T,i are summed over all triangles of the cell and transferred to the nodes in the direction perpendicular to the opposite vertex edge of that node leading to increase (or decrease) the area of one triangle at the expense of a decrease (or increase) of its neighbor triangles. We chose k A such that the influence on the cortex elasticity is minimal, yet triangle deformations are minimized. Note, these forces do not induce shear elastic effects in the network. Finally, the interactions between neighboring cells are accounted for by introducing repulsive forces ( F rep,i ) and the adhesive forces ( F adh,i ) between nodes belonging to different cells. The interaction forces were obtained as detailed in [4,20]. The model allows simulating the interaction between two arbitrarily shaped triangulated bodies but also between a triangulated body and non-triangulated simple geometric objects as spheres and flat substrates.
Importantly, the parameters of the spring network in Equation 0.3 can be related to macroscopic elastic constants, approximating the cell cortex by a thin shell. For the six fold symmetric triangulated lattice on the cell surface, the linear spring constant k s can be computed from the Young modulus E cor of the cortex with thickness h cor by [14,16,18]: Similarly, the bending stiffness of the cortex can be approximated by where ν is the Poisson ratio (= 0.3 for an equilateral 2D network of linear springs).
"Global" calibration approach, only valid for confined spheroid The relationship between forces, inter-cellular volume fraction and pressure established in section "Local" calibration approach, needed for experiment II can also be computed following a simpler, global approach in case of a confined volume as in experiment I. If compression rates are sufficiently slow as is the case in that experiment, the cells in the spheroid can reorganize, distributing the stress isotropically and homogeneously over the cells in the spheroid. For the intercellular spaces in the capsule the global intercellular volume fraction int and the cellular fraction cells are related by int + cells = 1 with cells = V cells /V caps . We can now use cells to re-parameterize Equation 14 assuming cells are homogeneously distributed over the spheroid during compression (compare Figure 8A). The apparent contact stiffnessẼ i in that case increases for all cells equally with cells being the equivalent of 1 −d ij in the previous approach.
With both different calibration methods, local and global, the CBM simulation results closely follow the DCM curve ( Figure 8D). However, the curve using the global calibration method is smoother than for the local calibration method as it represents an average over all cells thus disregarding local fluctuations. For high pressures (p > 2 kPa), both curves become nearly parallel. On the other hand, as the local approach does not require the existence of an enclosed (capsule) volume it can be used more generally as a cell-cell interaction force upon volumetric compression in many configurations as in experiment II.
We finally point out that no explicit representation of ECM has been considered in the model based on the observations the ECM fraction is rather small (personal communication) and approximately homogeneously and isotropically distributed in the intercellular spaces.

Model parameter and algorithm sensitivity
We also studied the potential influence of parameters calibrated for the free growth conditions which could not be directly inferred from the experiments. For the friction terms, the cell-cell friction, cell -ECM friction and cell-capsule friction were varied between 1 % ("low") and 200 % ("high") of their reference values. The effect is shown in Figure BA and indicates that even for low friction coefficients the results remain largely unaffected. For the cell-cell adhesion energy W and cell motility coefficients D, which were varied between 10 % and 500 % [6] , we did not observe any significant changes ( Figure BE-F). The strength of cell-cell adhesion has been shown to play a role in detachment, but to be of minor importance for multicellular systems under compression (e.g. [11]). Note, that also in MCSs growing in absence of externally applied forces cells are moderately compressed [10].
In conceptual analogy to experimental statistical procedures, we have performed growth expansion simulations in the thin capsule with the optimal parameter set but four different random seeds for the cell-specific cycle time, Young modulus and necrotic pressure threshold to test the effect of stochasticity (see Table 2, Figure BB). Even after more than 7 days of simulation, only very little mutual differences can be observed, while the slopes of the curves are the same. This can be attributed to self-averaging effects such that the variations on the level of individual cells cancel out at the population level (e.g. ref. [7]).

Comparison of calibration methods
For the local calibration procedure, we used the following constants in Equation 14 assuming strain hardening: a 0 = 0.4454 · 10 5 , a 1 = 2.347 · 10 5 , a 2 = −7.918 · 10 6 , a 3 = 6.615 · 10 8 , a 4 = −1.206 · 10 9 , a 5 = −3.091 · 10 9 , a 6 = 1.1239 · 10 10 , yielding a force cell -intercellular distance curve that matches well with the one obtained from the deformable cell simulations experiment (Figure 8). The function forẼ(d ij ) is monotonically increasing ( Figure BD). For the global calibration approach, we obtained a similar curve. In test simulations for the MCS growth in the thin capsule using contact inhibited growth for both local and global calibration, we found only a mutual variation of 5 − 7%, see Figure BC. On the other hand, if no calibration was assumed, a largely unrealistic capsule dilation was obtained ( Figure BC), and a consistent relation between pressure and cell density could not be established.

Influence of cell division algorithms
As far as either the volume or the duration of the cell in the cell cycle have passed a threshold value, we replace the mother cell by two daughter cells which are placed very close to each other [17,22]. This algorithm is different from the approach where cells slowly deform into dumbbells in mitosis phase before splitting into two daughter cells [11,21]. In our model two daughter cells are created next to each other instantaneously. During the mitosis period (which takes about 1.5h) the daughter cells do not grow. However, the two newly created cells can generate short-time artificial pressure peaks which disappear during the division time course due to small local re-arrangements. To test the impact of the pressure peak formation on our growth simulation results we implemented an smoothing algorithm that reduces these peaks by ensuring that (1) the mother cell divides in the direction of the least stress as derived from the local stress tensor (see section Measures for stress and pressure), and (2) a local energy minimum is sought by varying the distance between the daughter cells and computing the interactions with the other cells. While the smoothing algorithm reduces the short-time pressure peaks, we did not see significant differences in the results compared to simulations where this algorithm had been dropped.  Figure 3 (which fails to explain the data in experiment II, see Figure 7), and a linear function ("Linear -II") optimized to match the CT26 spheroid expansion in Experiment II (but here fails to reproduce the thin capsule data). (H) Data plots of MCS growth in a thick large capsule (∼ 400 µm) growth. The stress-dominated growth regime is too short to identify the stress-response.

5/9
Future extension: possible ways to include extracellular matrix An important future extension would be the inclusion of extracellular matrix (ECM), which is why we briefly outline how this can be implemented within the current model framework.
There are three approaches with different degree of detail and effort of implementation, one requiring an explicit model of ECM, while two of them might be readily implemented with modest effort.
(i) Approach 1: A sophisticated way would be to model the ECM explicitly e.g., as a network of visco-elastic springs (e.g. Ban et al., 2018 [23]). Such an approach is expensive both implementationwise and in terms of computing power. Alternatively, the effect of ECM can be estimated within two approximations within the framework of the established model.
(ii) Approach 2: One way is to perform the same simulation and calibration but now assuming that there is an ECM substance present in the capsule. To compute the influence of ECM of compression modulus K ecm , one could assume that the ECM is compressed homogeneously and isotropically (shear elastic effects with cells are neglected). In a next step the volume occupied by the (DCM) cells and ECM within the capsule would be estimated prior to compression of the capsule simulation (see Figure 8A) to obtain a corrected relation between cell-cell interaction force and cell-cell distance ( Figure 8C) that includes the compression and deformation forces of the ECM. For the Dextran experiment, the volume of the MCS could be estimated by the volume inside an elastic membrane enclosing the MCS and touching its surface.
The pressure necessary for ECM compression could be calculated from Equation 8 whereby dV then denotes the deviation of the ECM volume from its uncompressed volume as a consequence of compression leading to a relation p(V ) for the ECM alone. In the compression simulation used to determine the intercellular forces with cell-cell distance, the reduction of cell-cell distance would now cause both a compression of ECM by reducing its volume, and a deformation of the cells. Both, ECM and cell now resist the compression whereas in neglecting ECM, only the resistance of the cells was traced. As the volume of ECM can be calculated at any time from the difference of capsule volume and total volume of all cells, the pressure generated by compression of ECM during the compression simulations can be calculated, and taken into account in the calculation of the force on the cell by adding a force to each node that is proportional to the area to the node's associated cell surface area. This extends the curves in Figure 8C,D by the effect of ECM. For small compression modulus the force would be expected to be small. If at the same time, the compression modulus of the cell would be much larger than order (kPa) (see the discussion of cell compression moduli in the introduction), the repulsive branch of the curves in Figure 8 would become steeper. Approach 2 is a global approximation as it needs information from the entire MCS (the MCS or capsule volume) to estimate the effect of ECM and the simplest approximation.
(iii) Approach 3: The approximation could be based on a composite material approach in which the object that in the main body of the paper mimics the deformable cell, would instead be associated with a tissue unit composed of the cell plus a shell of ECM (for the concept in agent-based models, see Drasdo et. al. (2007) [11]). To determine the relation between force and cell-cell distance in such an approach, the radius of the object now composed of cell and ECM shell would be chosen larger than the cell radius alone to reflect the observed volume percentages of cells and ECM inside the MCS. The Young modulus of the cell cortex in the current approach would be replaced by an effective Young modulus taking into account the Young modulus of the ECM and the Young modulus of the cell cortex. As a simple illustrative example in point mechanics, it is conceptually equivalent to replacing two linear springs of spring constants k 1 (e.g. for the cell cortex) and k 2 (e.g. for the ECM) connected in series by one effective spring with spring constant k = 1/(1/k 1 + 1/k 2 ) (e.g. for cell cortex plus ECM). Similarly, the compression modulus of the cortex (that emerges from its Young modulus and Poisson ratio) can be replaced by an effective compression modulus, that other than the spring constants also contain the thicknesses of ECM and the cellular cortex. The effective compression modulus again contains information on compression moduli and volumes. This makes it possible to use the same approach as in the body of this paper by only re-parameterization of the cortex variables and the cell radius. However, in a configuration where cells touch each other, the ECM shell is already under slight compressive stress, which in reality might not be the case. From the definition of the composite material parameters and the geometrical quantities it is possible to track the position of the cell surface (such as its is possible to calculate in the above point mechanics example the contact point between the two springs). In the extremal case where the cell (including cortex) compression and Young moduli would be much larger than the compression modulus of the ECM, deformation and compression of a cell / ECM unit would basically extend solely on the ECM shell. Approach 3 would be local and is expected to give reasonable approximations in particular if the compressibility of the ECM is much bigger than that of a cell. Performing simulations with a DCM with composite parameters would lead to a different relationship of interaction force versus distance, which then could again be used within CBM simulations.