Exact Solutions of Coupled Multispecies Linear Reaction–Diffusion Equations on a Uniformly Growing Domain

Embryonic development involves diffusion and proliferation of cells, as well as diffusion and reaction of molecules, within growing tissues. Mathematical models of these processes often involve reaction–diffusion equations on growing domains that have been primarily studied using approximate numerical solutions. Recently, we have shown how to obtain an exact solution to a single, uncoupled, linear reaction–diffusion equation on a growing domain, 0 < x < L(t), where L(t) is the domain length. The present work is an extension of our previous study, and we illustrate how to solve a system of coupled reaction–diffusion equations on a growing domain. This system of equations can be used to study the spatial and temporal distributions of different generations of cells within a population that diffuses and proliferates within a growing tissue. The exact solution is obtained by applying an uncoupling transformation, and the uncoupled equations are solved separately before applying the inverse uncoupling transformation to give the coupled solution. We present several example calculations to illustrate different types of behaviour. The first example calculation corresponds to a situation where the initially–confined population diffuses sufficiently slowly that it is unable to reach the moving boundary at x = L(t). In contrast, the second example calculation corresponds to a situation where the initially–confined population is able to overcome the domain growth and reach the moving boundary at x = L(t). In its basic format, the uncoupling transformation at first appears to be restricted to deal only with the case where each generation of cells has a distinct proliferation rate. However, we also demonstrate how the uncoupling transformation can be used when each generation has the same proliferation rate by evaluating the exact solutions as an appropriate limit.


Numerical Methods
We will now describe how we obtain numerical solutions to the system of coupled equations given by on 0 < x < L(t), with zero diffusive flux boundary conditions, at both x = 0 and x = L(t), for each species. The first step is to apply the boundary fixing transformation ξ = x/L(t), and recall that v = ξdL(t)/dt, which allows us to re-write the system as ∂C 4 ∂t = D L 2 (t) on 0 < ξ < 1, with ∂C i /∂ξ = 0 at ξ = 0 and ξ = 1 for i = 1, 2, 3 and 4. We discretize Equations (5)-(8), with a central finite difference approximation, on 0 < ξ < 1, with uniform mesh spacing δξ. The resulting system of coupled ordinary differential equations is integrated through time using a backward Euler approximation with uniform time steps of duration δt. The numerical solutions are integrated through time sequentially, treating Equations (5)-(8) separately. The resulting system of tridiagonal linear equations for each component is solved using the Thomas algorithm.
All numerical results in the main manuscript are obtained for choices of δξ and δt that produce grid-independent results. The numerical solutions on the fixed domain, 0 < ξ < 1, are then re-expressed on the growing domain by applying the inverse boundary fixing transformation, x = ξL(t).

2/11
Combined homogeneous Dirichlet and homogeneous flux boundary conditions In the main part of the manuscript we demonstrate how to solve the governing equations for homogeneous Neumann (no-flux) boundary conditions. It is also possible to solve problems that involve a combination of homogeneous Dirichlet and zero flux boundary conditions, and we will demonstrate this by considering on 0 < x < L(t) with ∂a/∂x = 0 at x = 0, and a = 0 at x = L(t). This means that we have zero diffusive flux at the origin, and a homogeneous Dirichlet condition at the moving boundary x = L(t). We proceed in the usual way by applying the boundary fixing transformation, ξ = x/L(t), and then re-scaling the time variable where λ n = (2n − 1)π/2. In this form, the Fourier coefficients, Ψ n , can be chosen so that the solution matches the initial condition for a(ξ, 0). Once we determined the Fourier coefficients we can apply the inverse transforms to give a(x, t), as required.
explain how this procedure can be adapted to deal with homogeneous and nonhomogeneous Dirichlet boundary conditions. However, the adaption of our technique to deal with nonhomogeneous flux boundary conditions is not standard, and we explain the key steps here. If we are interested in solving a coupled system of reaction-diffusion equations on a growing domain with nonhomogeneous flux boundary conditions, we first apply the Sun and Clement transformation to uncouple the system. This will lead to a systems of uncoupled partial differential equations.
Each equation in the system has the form where where choosing satisfies the nonhomogeneous boundary conditions for a(ξ, T ). After making these with homogeneous boundary conditions ∂v/∂ξ = 0 at ξ = 0 and ∂v/∂ξ = 0 at ξ = 1.
Here, h(ξ, T ) is known, and given by To proceed we assume that both the known function, h(ξ, T ), and the unknown function, v(ξ, T ), can be written in terms of their cosine eigenfunction expansions where To solve for v 0 (T ) and v n (T ), we substitute Equations (21)-(22) into Equation (17), which gives an infinite system of ordinary differential equations whose solution can be obtained by using an integrating factor For exponential growth we have σ(t) = α and To simplify the algebra we define It is important to note that as t → ∞ we have T (t) → (1/κ) − , so that 1 − κT > 0 for all t > 0. This means that f (T ), g(T ) and w(ξ, T ) are all well-defined for all t > 0.
With the simplest possible initial condition, a(x, 0) = 0, we obtain v 0 (T ) = 2Ω 3(κ + 2ν) where E s (z) = ∞ 1 e −zt /t s dt is the exponential integral [1]. With this information the solution can be written as which can be rewritten in the original coordinate system, a(x, t), by recalling that T (t) = (1 − e 2αt D)/(2αL 2 (0)) and ξ = x/(L(0)e αt ). Although the eigenfunction expansion technique described here produces an exact solution for an exponentially growing domain with a simple initial condition, a(x, 0) = 0, we have found that applying the technique to other choices of L(t) and a(x, 0) can lead to intractable integrals in the definition of v 0 (T ) and v n (T ).

8/11
Approximating nonlinear models by linearised models A common theme in the Introduction and Discussion sections of the main document is that one of the key limitations of our work is that the exact solution strategy is valid only for linear partial differential equations whereas many reaction-diffusion models used to represent collective cell spreading processes are nonlinear. We will now demonstrate the ability of a linear reaction-diffusion equation to approximate the solution of a related nonlinear reaction-diffusion equation on a growing domain. To demonstrate this we will consider a single, uncoupled, reaction-diffusion equation on a growing domain, where D is the diffusivity, v is the advective velocity associated with domain growth, k is the proliferation rate and S(C) is a generalised source term. Setting S(C) = C means that we are considering a linear source term whereas setting S(C) = C(1 − C) is a logistic-type nonlinear source term. This kind of nonlinear source term is often used to mimic crowding effects since we have S(1) = 0 which means that the proliferation rate becomes zero when the density increases to the carrying capacity density, which here we have set to be unity. To demonstrate the differences between setting S(C) = C and S(C) = C(1 − C) we will obtain numerical solutions of Equation (35). When we set S(C) = C(1 − C), numerical solutions of the nonlinear equation are obtained by using an iterative fixed-point method [3].