Mechanical Models of Pattern and Form in Biological Tissues: The Role of Stress–Strain Constitutive Equations

Mechanical and mechanochemical models of pattern formation in biological tissues have been used to study a variety of biomedical systems, particularly in developmental biology, and describe the physical interactions between cells and their local surroundings. These models in their original form consist of a balance equation for the cell density, a balance equation for the density of the extracellular matrix (ECM), and a force-balance equation describing the mechanical equilibrium of the cell-ECM system. Under the assumption that the cell-ECM system can be regarded as an isotropic linear viscoelastic material, the force-balance equation is often defined using the Kelvin–Voigt model of linear viscoelasticity to represent the stress–strain relation of the ECM. However, due to the multifaceted bio-physical nature of the ECM constituents, there are rheological aspects that cannot be effectively captured by this model and, therefore, depending on the pattern formation process and the type of biological tissue considered, other constitutive models of linear viscoelasticity may be better suited. In this paper, we systematically assess the pattern formation potential of different stress–strain constitutive equations for the ECM within a mechanical model of pattern formation in biological tissues. The results obtained through linear stability analysis and the dispersion relations derived therefrom support the idea that fluid-like constitutive models, such as the Maxwell model and the Jeffrey model, have a pattern formation potential much higher than solid-like models, such as the Kelvin–Voigt model and the standard linear solid model. This is confirmed by the results of numerical simulations, which demonstrate that, all else being equal, spatial patterns emerge in the case where the Maxwell model is used to represent the stress–strain relation of the ECM, while no patterns are observed when the Kelvin–Voigt model is employed. Our findings suggest that further empirical work is required to acquire detailed quantitative information on the mechanical properties of components of the ECM in different biological tissues in order to furnish mechanical and mechanochemical models of pattern formation with stress–strain constitutive equations for the ECM that provide a more faithful representation of the underlying tissue rheology. Supplementary Information The online version contains supplementary material available at 10.1007/s11538-021-00912-5.

In this document (supplementary information) we report on details of the numerical schemes employed to obtain the numerical results presented in the above named publication (main publication). For details of the equations referred to in this document, please see the main publication. Numerical solutions are computed using Matlab, and the files containing the code corresponding to the schemes presented below can be found on GitLab (https://git-ce.rwth-aachen.de/alf.gerisch/VillaEtAl2021BullMathBiol).
Numerical schemes for the system of PDEs (10), (11) and (16) Numerical solutions for the system of implicit, time-dependent and spatially one-dimensional PDEs (10), (11) and (16) are obtained exploiting the Method of Lines. We make use of a uniform discretisation of the spatial domain [l, L] consisting of K + 1 grid points, or grid cell centres, while, at first, leaving the time variable continuous. We denote the spatial grid width by ∆x. The normalised cell density n(t, x), the normalised ECM density ρ(t, x) and the displacement of a material point of the cell-ECM system u(t, x) are approximated as Thanks to the periodic boundary conditions we have N 0 (t) = N K (t) , P 0 (t) = P K (t) and U 0 (t) = U K (t) , and consequently have 3 × K time-continuous approximations to determine. We collect them in the vectors N (t), P (t), U (t) and denote their time-derivatives by N (t), P (t), U (t). The discretization of the spatial derivatives in the PDE system will then result in an implicit system of 3 × K ODEs for the variables N (t), P (t), U (t) and their time-derivatives of the following form (S.1) In this system f n (N, P, N , U ) = 0, f ρ (P, P , U ) = 0 and f u (N, P, U, N , P , U ) = 0 are each systems of K ODEs obtained, respectively, from PDEs (10), (11) and (16), using second-order central finite difference approximations for the spatial derivatives and the first-order upwind scheme for the advection terms, as detailed below for each equation.
In order to solve system (S.1), we make use of the Matlab solver ode15i, which uses a variable-order (orders 1 to 5) backward difference formula (BDF) method in a form suitable to an implicit system of ODEs. Initial conditions N (0), P (0) and U (0) are given by the appropriate equivalent of initial conditions (28), and we make use of the Matlab function decic to obtain consistent initial conditions N (0), P (0) and U (0) such that (S.1) is satisfied at initial time t = 0.
Useful matrices. In order to apply the first-order upwind scheme we need to compute variables and derivatives at the grid cell interfaces, i.e. half-way between grid points, in addition to those at the grid cell centres. We here clarify the notation adopted throughout the rest of this document. The K×K matrices M x and M xx are used to approximate, using second-order finite differences, the first-order and the second-order derivatives in space, respectively, of a periodic grid function at the grid cell centres and are therefore given by We make use of the notation − · to indicate a shift from the grid cell centres to the (right) grid cell interfaces. In particular to approximate the value of a periodic grid function at these grid cell interfaces we multiply it by the K×K matrix In addition, the K×K matrices − M x and − M x are used to approximate the first-order derivatives in space of a periodic grid function at the (right) grid cell interfaces, when the grid function is given in the grid cell centres, and at the grid cell centres, when the grid function is given in the (right) grid cell interfaces, respectively. These are given by Note that, even though these two matrices are multiplied by 1/∆x, they still stem from second-order finite difference approximations, calculated on a staggered grid shifted by half the grid cell width. Convention: In the formulas which follow below, we use the convention that any product of a matrix from above with a vector of length K is a matrix-vector product but any operation between two vectors, in particular multiplication, division, or exponentiation, are understood element-wise.
Numerical scheme for the balance equation (10). We rewrite the balance equation (10) as which, upon spatial discretisation, leads to the following system of K ODEs where the matrix − M x is defined in (S.4) and the advective flux − F at the grid cell interfaces is computed using first-order upwinding, i.e. Numerical scheme for the transport equation (11). We rewrite the transport equation (11) as which, upon spatial discretisation, leads to the following system of K ODEs Numerical scheme for the force-balance equation (16). We solve the system of PDEs (10), (11) and (16) for the Kelvin-Voigt (3) and the Maxwell (4) models. In these cases we have b 2 = a 2 = 0, and the force-balance equation (16) reads as Upon spatial discretisation, this leads to the following system of K ODEs where I is the K×K identity matrix and M x is defined in (S.2). This scheme is valid as long as b 2 = a 2 = 0 and can therefore also be applied when considering the linear elastic model (1), the linear viscous model (2), and the SLS model (5). On the other hand, in the case where b 2 = 0 (i.e. when the Jeffrey model (6) is considered) the above numerical scheme cannot be directly employed due to the presence of a second-order derivative in t. We could still, however, take a similar approach and make use of the ode15i solver by introducing extra variables for the first-order derivatives in t of n and ρ, thus formally reducing the PDE (16) to first-order in time, at the cost of increasing the number of equations in the Method of Lines ODE system.

Numerical scheme for the system of PDEs (29)
Similarly as done for the spatially one-dimensional model, numerical solutions for the system of implicit, time-dependent and spatially two-dimensional PDEs (29), together with (30)-(32), are obtained exploiting the Method of Lines. We make use of a uniform discretisation of the square spatial domain [l, L] × [l, L] consisting of (K + 1)×(K + 1) grid points, while leaving the time variable continuous. The spatial grid width, in both spatial directions, is denoted by ∆x again. The normalised cell density n(t, x 1 , x 2 ), the normalised ECM density ρ(t, x 1 , x 2 ) and the displacement of a material point of the cell-ECM system u(t, x 1 , x 2 ) = u 1 (t, x 1 , x 2 ), u 2 (t, x 1 , x 2 ) are approximated as Thanks to the periodic boundary conditions, we can drop the index values i = 0 and j = 0 and consequently have 4×K 2 time-continuous approximations to determine. We collect them in the matrices N (t), P (t), U 1 (t), U 2 (t) and denote their time-derivatives by N (t), P (t), U 1 (t), U 2 (t). The discretization of the spatial derivatives in the PDE system will then result in an implicit system of 4 × K 2 ODEs for the variables N (t), P (t), U 1 (t), U 2 (t) and their time-derivatives of the following form In this system f n (N, P, N , U 1 , U 2 ) = 0, f ρ (N, P, N , U 1 , U 2 ) = 0, f u1 (N, P, U 1 , U 2 , N , P , U 1 , U 2 ) = 0 and f u2 (N, P, U 1 , U 2 , N , P , U 1 , U 2 ) = 0 are each systems of K 2 ODEs obtained from the system of PDEs (29), using second-order central finite difference approximations for the spatial derivatives and the first-order upwind scheme for the advection terms, as detailed below for each equation.
In order to solve system (S.15), we make, similarly to the spatially one-dimensional case, use of the Matlab solver ode15i. Initial conditions N (0), P (0), U 1 (0) and U 2 (0) are given by the appropriate equivalent of initial conditions (36), and we make use of the Matlab function decic to obtain consistent initial conditions N (0), P (0), U 1 (0) and U 2 (0) such that (S.15) is satisfied at initial time t = 0.
Useful functions In order to solve the system (S.15) we need to compute variables and derivatives at the grid cell centers and interfaces, both in the x 1 -and the x 2 -direction. We here introduce the functions that will be used in the rest of this document to compute the aforementioned quantities in the different directions. These rely on the fact that the matrices (S.2)-(S.4) act on column vectors and therefore, when applied to an K×K argument matrix, they will act on each column of that, which in our framework corresponds to computing the quantity of interest in the x 1 -direction. In order to compute the same quantities in the x 2 -direction, we need the operating matrix to act on each row of the argument matrix of interest, which can be achieved by matrix transposition of the argument matrix before and of the product matrix after matrix multiplication. Hence the functions M x1 (N ) and M x2 (N ) are used to approximate the first-order derivative of the variable of interest, say N , at the grid cell centres in the x 1and x 2 -directions respectively, and are defined as with the function A − v 1 , N given by (S.7) together with definitions (S.8) and (S.9). Convention: With the definitions above, we have hidden all applications of the matrices from the spatially one-dimensional case in newly defined functions. Consequently, in the formulas which follow below, we use the convention that any further operation between matrices, in particular multiplication, division, or exponentiation, are understood element-wise.