A Continuum Mathematical Model of Substrate-Mediated Tissue Growth

We consider a continuum mathematical model of biological tissue formation inspired by recent experiments describing thin tissue growth in 3D-printed bioscaffolds. The continuum model, which we call the substrate model, involves a partial differential equation describing the density of tissue, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{u}}(\hat{{\mathbf {x}}},{\hat{t}})$$\end{document}u^(x^,t^) that is coupled to the concentration of an immobile extracellular substrate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{s}}(\hat{{\mathbf {x}}},{\hat{t}})$$\end{document}s^(x^,t^). Cell migration is modelled with a nonlinear diffusion term, where the diffusive flux is proportional to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{s}}$$\end{document}s^, while a logistic growth term models cell proliferation. The extracellular substrate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{s}}$$\end{document}s^ is produced by cells and undergoes linear decay. Preliminary numerical simulations show that this mathematical model is able to recapitulate key features of recent tissue growth experiments, including the formation of sharp fronts. To provide a deeper understanding of the model we analyse travelling wave solutions of the substrate model, showing that the model supports both sharp-fronted travelling wave solutions that move with a minimum wave speed, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c = c_{\mathrm{min}}$$\end{document}c=cmin, as well as smooth-fronted travelling wave solutions that move with a faster travelling wave speed, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c > c_{\mathrm{min}}$$\end{document}c>cmin. We provide a geometric interpretation that explains the difference between smooth and sharp-fronted travelling wave solutions that is based on a slow manifold reduction of the desingularised three-dimensional phase space. In addition, we also develop and test a series of useful approximations that describe the shape of the travelling wave solutions in various limits. These approximations apply to both the sharp-fronted and smooth-fronted travelling wave solutions. Software to implement all calculations is available at GitHub. Supplementary Information The online version contains supplementary material available at 10.1007/s11538-022-01005-7.

S1 Numerical methods 1 S1.1 Two-Dimensional partial differential equations 2 Numerical solutions of the substrate model in two-dimensions are obtained by re-writing Equations (1)-(2) as where we have written the mathematical model in terms of a general nonlinear 3 diffusivity function D(s), and general source terms, f (u) and g(u, s). To be con-4 sistent with our dimensional model we have D(s) = Ds/K s , f (u) = λu(1 − u/K u ) 5 and g(u, s) = r 1 u − r 2 s, but our numerical method can deal with other functional 6 forms if required. Our aim is to obtain numerical solutions of Equations (S1)-(S2) 7 on the square domain Ω = {(x, y), 0 < x < L, 0 < y < L}. For convenience we 8 assume that the origin is at the lower left corner of the domain and we discretise 9 Ω on a spatially uniform finite difference mesh with mesh spacing ∆x = ∆y > 0. 10 We index the mesh in the usual way so that the coordinates of each mesh point 11 are (x i , y j ), with i = 0, 1, 2, . . . , I and j = 0, 1, 2, . . . , J. Since we always consider a 12 square mesh we have I = J. All numerical results correspond to a 101 × 101 mesh 13 which, with L = 300 µm, gives ∆x = 3 µm. We found that solutions obtained on 14 a finer mesh gave visually indistinguishable results for the parameter values that 15 we considered.
We solve Equations (S1)-(S2) using a standard method of lines approach so that at each internal mesh point we have where we have approximated the internode diffusivity with an arithmetic average.

17
These discretised equations are valid at central nodes, i = 1, 2, . . . , I − 1 and To solve Equations (3)-(4) we consider a domain 0 < x < L that we discretise into m equally-sized intervals with spacing ∆x. We approximate Equations (3)-(4) using a central difference approximation for the spatial derivative. Since we use this algorithm to study long-time travelling wave solutions we approximate the temporal derivative with an implicit Euler approximation, giving We estimate the travelling wave speed by specifying a particular contour value, 37 u(x, t) = u * and use linear interpolation to estimate x * such that u(x * , t) = u * at 38 each time step. With this data we then calculate at each time step, which we find settles to a constant value for sufficiently large t.

40
Given this time series of estimates for c we fit a straight line to the late-time data to provide an estimate of c. All results in this work correspond to u * = 0.5, but we 42 find that our results are insensitive to this choice and other values of u * ∈ (0, 1) 43 give the same results provided that δt and h are chosen to be sufficiently small. S1.4 Phase plane on the slow manifold 45 We Here we provide numerical evidence to test the hypothesis that the shape of 52 smooth-fronted travelling waves are given by Equations (19)-(21). Results in Fig-53 ure S1(a) show a smooth-fronted travelling wave with c = 2.00 for r 1 = r 2 = 1.

54
Using our long-time numerical PDE solution we plot U (z), S(z) and W (z) at 55 the leading edge of the travelling wave in Figure S1(b). The profiles in Figure   56 S1(b) are a magnified view of the region contained within the purple rectangle in 57 Figure S1(a). At the scale shown in Figure S1(b) we clearly see U (z), S(z) and 58 W (z) decaying to zero with z, and each numerical profile is superimposed with For (e)-(f) the initial conditions are given by Equations (6)-(7). All numerical PDE solutions correspond to ∆x = 1 × 10 −4 , ∆t = 1 × 10 −3 and = 1 × 10 −4 . S3 Desingularised phase space and slow manifold reduction. All results correspond to R = 2 (r 1 = 1, r 2 = 0.5). Results in: (a)-(b) correspond to a smooth-fronted travelling wave with c 2 = 10; (c)-(d) correspond to a smooth-fronted travelling wave with c 1 = 1; and, (e)-(f) correspond to a sharp-fronted travelling wave with c min = 0.30. Results in the left-most column show the three-dimensional desingularised phase space with the invaded equilibrium point (green dot), the uninvaded equilibrium point (blue dot) and the slow manifold (grey surface). Results in the right-most column show the vector field on the slow manifold, superimposed with several solution trajectories, including the heteroclinic orbit (blue) and several unphysical trajectories (red). The heteroclinic orbit is obtained by solving Equations (3)-(4) numerically with appropriate initial conditions. For (a)-(b) and (c)-(d) the initial conditions are given by Equations (8)-(9) with a = 1/10 and a = 1, respectively. For (e)-(f) the initial conditions are given by Equations (6)-(7). All numerical PDE solutions correspond to ∆x = 1 × 10 −4 , ∆t = 1 × 10 −3 and = 1 × 10 −4 .

S2.4 Phase space for c < c min 74
In this section we explore the consequences of setting c < c min in the phase space.

75
Results in Figure S4   In Equation (52) we left the expression for U 1 (ẑ) as an integral. Here, we give present a solution for U 1 (ẑ) for some special choices. For r 2 = 1 we obtain and when r 2 = 2 we obtain where Li 2 (x) is a special function called the dilogarithm function that is given