Free boundary problems for Stokes flow, with applications to the growth of biological tissues

We formulate, analyse and numerically simulate what are arguably the two simplest Stokes-ﬂow free boundary problems relevant to tissue growth, extending the classical Stokes free boundary problem by incorporating (i) a volumetric source (the nutrient-rich case) and (ii) a volumetric sink, a surface source and surface compression (the nutrient-poor case). Both two- and three-dimensional cases are considered. A number of phenomena are identiﬁed and characterised thereby, most notably a buckling-associated instability in case (ii).


Mathematical approaches to modelling tissue growth
Macroscale modelling of the growth of biological tissue raises a wide variety of issues in the formulation and mathematical treatment of free boundary problems as well as providing insight into questions of practical significance. The earliest work in this direction focussed on cancer growth (see [13], for example, and [4] for a more recent overview) while recent developments in tissue engineering provide further motivation (e.g., [22]), as do questions in developmental biology (e.g., [2,16,24]).
Key phenomena in tissue growth are associated with biomechanics (notably, the stresses induced by cell division within the tissue) as well as with nutrient transport and consumption (only those cells that are adequately supplied undergo cell division, while those that are under-resourced are quiescent or subject to cell death). Here, we focus explicitly on some of the simplest mathematical descriptions of both of these aspects, seeking in particular to clarify the mechanisms that, whilst present in models that aspire to greater biological realism, are at times obscured by the levels of complexity. In terms of mechanics, following [13], a Darcy constitutive relationship has very often been adopted; while appropriate for tissue-engineering applications of growth in a porous scaffold, this seems inappropriate both in the context of in vivo applications and for many in vitro ones such as the (unconstrained) growth of tumour spheroids or embryoid bodies where there is no static matrix through which the tissue is flowing. A recent study, [9], follows a similar approach to the current paper in exploring some of the properties of the simplest Darcy-based models, but here we instead adopt, for the reasons highlighted above, a Stokes-flow constitutive assumption, as in [10], [19] and [21]. We remark that extensive current work investigates elastic effects (see for example, the book by [12]), which raise significant additional challenges -such as those associated with assigning a reference configuration to tissue newly generated by cell division. We eschew such complications here in keeping with our intention to explore models of minimal complexity that never the less retain essential features. Thus, we are in effect assuming that stress relaxation occurs over timescales significantly shorter than those of growth.

Structure of the manuscript
The remainder of the manuscript is structured as follows; in 1.3 we discuss some numerical methods for the approximation of free boundary problems in fluid flows. In 2 we derive a onephase Stokes-flow free boundary problem from a simple model for tissue growth considering two distinct settings, namely nutrient-rich and nutrient-poor cases, and we discuss some special cases in which analytical progress is possible. In 3 we present a weak formulation of the free boundary problem and a finite element approximation based on this weak formulation. In 4.1 we present simulations of the Stokes-flow free boundary problem in the nutrient-rich case in 2d and 3d . In 4.2 we present simulations of the Stokes-flow free boundary problem in the nutrient-poor case in 2d and 3d . Finally, in 5 we summarise our results and discuss an avenue for further investigation.

Numerical methods for Stokes-flow free boundary problems
The structure of the model for tissue growth that we derive corresponds to a one-phase Stokes-flow free boundary problem with surface tension. A large number of approaches have been developed for the approximation of free boundary problems appearing in fluid flows ( [15]), often accounting for multiple phases and considering the full Navier-Stokes system. Broadly speaking, the methods may be split into two classes, implicit methods in which fluid interfaces are not approximated directly, such as the volume of fluid (VOF) methods ( [17]) in which an indicator function is introduced to track the region occupied by fluid phases (see [11] for a review), level set methods in which a function is introduced whose zero level-set corresponds to the interface between phases (e.g., [25]) and diffuse interface or phase-field methods that introduce a phase-field function to track the evolving phases, the interface between phases being approximated via a diffuse region (e.g., [18]). The aforementioned approaches all have the advantage of being robust to large deformations and they allow for topological changes in the interface. However, the accurate resolution of the interfacial regions requires a fine grid and this can lead to challenges in terms of computational cost. The other class of methods consists of parametric and front tracking methods which seek to approximate the interfacial regions directly by fitting a mesh to the interface and evolving this surface mesh along with the fluid (e.g., [23] and [1]). Such methods allow accurate computation of the interface, although it can be challenging to deal with large deformations due to degeneration of the surface mesh; moreover, topological changes of the interface are in general not permissible. For the present study, since the simple one-phase Stokes-flow free boundary problem we consider does not precipitate any of the issues that make the application of the parametric approach challenging, we employ a parametric finite element method for the approximation. We note that one extension of the present work we wish to consider is to allow for a more detailed model of the interfacial FREE BOUNDARY PROBLEMS FOR STOKES FLOW 435 dynamics which in itself is well suited to approximation via parametric methods such as the surface finite element method ( [8]). We remark that in the 2d case, complex variable theory can be exploited to develop highly efficient and accurate methods for the approximation of free boundary problems of the form considered in this work (e.g., [3]); however, we wish to adopt an approach that is not restricted to the 2d case.
2. Derivation of the Stokes-flow free boundary problem and investigation of some special cases 2.1 Asymptotic derivation of the one-phase Stokes free boundary problem We derive the model pursued subsequently from a simple two-phase framework in which the phases in question are cells (volume fraction n) and water (volume fraction w), with new cells being formed from water and water being released on death. The resulting mass continuity equations read 1 where S is the interconversion rate between phases (to be specified later -we note again that the simplest possible assumptions are being made here in terms of the composition of cells and of necrotic material) and u .n/ and u .w/ are the respective velocity fields of the two phases. Given that inertial effects are irrelevant, we take the momentum equations to be 0 D r n .n/ Á C p I rn where .n/ and .w/ are the intraphase stress tensors, p I is the interphase pressure and 1=k is the interphase drag coefficient (k being proportional to permeability). For the current purposes we can take the water to be inviscid, i.e.
.w/ ij D p w ı ij and treat the cell phase as a Newtonian fluid, so that where p n and p w are intraphase pressures and n and n are respectively the shear and bulk viscosities of the cell phase. It should be noted that the transfer S between phases induces a nonzero bulk viscosity term even though both phases are individually incompressible. If p n is identified with 3 tr .n/ , a choice often made, then relates n to n ; in any case the n term can for some purposes be absorbed into p n (cf. (2.6) below). To close the system, two expressions relating p n ; p w and p I are required, as well as boundary conditions (consideration of specific initial conditions will not be required). We shall, for simplicity, set p w D p I D P; p n D P C˙n.n/; this sufficing for the current purposes; the former introduces the notation P , while the extra pressure˙n embodies cell-cell interactions. On the interface, the simplest appropriate conditions are (adopting the summation convention over j and i in the second and third conditions) .n/ ij where 2 R denotes the surface tension of the interface, H denotes the mean curvature (sum of the principal curvatures) of the interface, with the convention that H > 0 if the interface is a spherical surface, D . 1 ; 2 ; 3 / T the outward-pointing normal to the interface and v the normal velocity of the interface.
Rescaling as appropriate to derive the one-phase model studied in the remainder of the paper, we rewrite˙n.n/ as˙n ..1 n/="/ and then set The rescaling of k implies easy transport of water through the tissue, enabling the problem for n to be decoupled, while the above assumption on˙n enables the rapid water flow required to sustain cell division even though w 1. Under these scalings, the leading order problem for the cell phase reads n D 1 together with r u .n/ D S; 0 D r .n/ ; .n/ with interface conditions as in (2.4), where p Á P C˙n.w/ n r u .n/ : (2.6) The system (2.5) is that studied below in two cases; S is typically made to depend on nutrient levels and we consider the two simplest possibilities -S a positive constant, corresponding to nutrient-rich conditions, and S a negative constant corresponding to a constant rate of cell death in a nutrientstarved necrotic core. In the former case, the above description applies as it stands; in the latter case, cell division is present in a thin rim (see, for example, the figures and data in [14,26]) so the situation is more complicated -an asymptotic analysis ( [21]) of the thin-rim limit, which collapses the rim onto the interface, shows that the conditions in (2.4) need to be modified to .n/ ij where 0 is a positive constant associated with compressive stresses arising from the cell division within the rim and q is a further positive constant, reflecting that in the limit cell division corresponds FREE BOUNDARY PROBLEMS FOR STOKES FLOW 437 to a surface source of material (this effect being analysed in the Darcy case in [9]). In (2.7) the case 0 > implies negative surface tension (net surface compression) with significant implications in what follows, wherein 0 6 will also be investigated. We recall that in the nutrient-rich case only positive surface tension is of interest, with q D 0: (2.5)-(2.7) allow p and u .n/ to be extracted given S and henceforth we discard the .n/ superscripts. Having solved for these, w can be obtained at leading order from with the boundary condition˙n .w/ D p C n S (2.9) giving w on the interface. (2.8)-(2.9) is a nonlinear elliptic problem for w with p determined from (2.5)-(2.7) and with k typically being another specified function of w. Because the interface evolution is independent of w, we shall not need to analyse this w formulation, though note thaṫ n .w/ needs to be such that w > 0 holds everywhere. The " ! 0 limit just described allowṡ n .w/ to pull in water rapidly in the nutrient-rich case to maintain growth with w 1, while for S < 0 water is transported to the rim from the core as well as from the exterior.
It is instructive to derive from (2.5) the most general p and u .n/ for which .n/ D 0 for uniform S, namely (in three dimensions) Sx C rigid body motions: (Radial growth, i.e., u .n/ / x, will reappear shortly.) For well-known reasons, the Stokes-flow prescription (2.5) can determine u .n/ only up to rigidbody motions; by reinstating inertia terms (which are otherwise entirely negligible) these could be determined in the usual way via solvability conditions on the overall linear and angular momenta of the tissue.

The Stokes-flow free boundary problem
The limiting problem derived in 2.1 forms the basis of our subsequent numerical investigations. As motivated above the problem of interest corresponds to the following free boundary problem where˝.t/ R d is the time dependent region consisting of the tissue, 2 R with > 0 denotes the viscosity of the tissue and 2 R can take either sign, according to whether a volumetric source or sink is under consideration. We denote by .t/ the time dependent interface that describes the moving boundary of the tissue and throughout we make the assumption that . .t// t2OE0;T is a smooth evolving hypersurface. We prescribe the initial condition .0/ D 0 : (2.11) The boundary conditions involve with Id 2 R d d denoting the identity matrix and the second equality defines the deformation tensor D.u/ D 1 2 ru C .ru/ T . We enforce the following boundary conditions on , denoting by and v the unit (outward) normal and normal velocity of where we recall that q denotes a surface source term, H denotes the mean curvature (sum of the principal curvatures) of with the convention that H > 0 if˝.t / is a sphere and 2 R the surface tension.
We can write the full free boundary problem (2.10)-(2.13) as (2.14) Following ( [7]), to approximate the mean curvature of , we make use of the identity where x denotes the identity map on .t/ and ∆ .t / denotes the Laplace-Beltrami operator on .t/ which is given by , the surface divergence of the surface gradient on .t/.

Special cases
We now discuss some simplified settings in which analytical progress can be made in terms of understanding the solutions to (2.14).

Radial symmetry.
Here we consider the one-dimensional, cylindrically symmetric and spherically symmetric growth of three-dimensional tissues. In each case the off diagonal components of the stress tensor are zero. Since we are taking to be constant, the results are as follows: (i) One-dimensional case (using Cartesian coordinates .x; y; z/ T , the free boundary corresponding to x D s.t /). For growth in the x direction, hold, with the other two components of the momentum equation being trivially satisfied. Hence, The pressure is thus uniform and, focussing on > 0, growth is prohibited in the y and z directions, with proliferation thereby inducing a compressive pressure (inclusion of a bulk viscosity term modifies its value in an obvious fashion), potentially inducing buckling. The free boundary satisfies P s D s (nutrient-rich); or P s D q C s (nutrient-poor): (ii) Cylindrically symmetric (with cylindrical polar coordinates . ; ; z/ T , the free boundary corresponding to D s.t/). In this case, Accordingly and the position of the free boundary is given by The proliferation pressure is lower than in (i) because the constraint on growth operates in the z direction only. (iii) Spherically symmetric (with spherical polar coordinates .r; Â; / T , the free boundary corresponding to r D s.t /). Here, Now, with the position of the free boundary given by Here there are no geometrical constraints on growth, surface tension/compression being the only source of non-zero stress components. Moreover, including bulk viscosity as in (2.3) leads to the same conclusion on p, whereas under the same assumption we have p D .4=3/ in case (i) and p D .1=3/ C =s in case (ii).
The solutions in the special cases (i), (ii) and (iii) can form the basis of linear-stability analysis (not reported here), a particular implication being that the < 0 case can be expected to be illposed, with arbitrarily short perturbations growing arbitrarily fast; the associated instability can be interpreted physically as a buckling instability ( [20]). The results of case (iii) can also be used to motivate the analysis in 2.3.2.
2.3.2 Zero surface tension. When D 0, significant progress is possible without imposing radial symmetry, provided unconstrained growth is considered: in this case the radial velocity component is given by with the other velocity components being zero, where d is the dimension both of the tissue and of the space in which it is growing (in contrast to the considerations in 2.3.1); since the stress components are all then zero, the boundary conditions (2.12) are automatically satisfied. Writing the free boundary as F .x; t/ D 0, in the nutrient-rich case the kinematic condition reads In consequence, ; i.e., the tissue grows without change of shape, as noted in ( [19]). In the nutrient-poor case, we instead have (we recall that < 0 holds in the nutrient-poor case) leads to F .X ; / being governed by @ F .X ; / D q jrF .X ; /j ; this being the interfacial-dynamics law for an interface evolving with constant outward-normal velocity q; in consequence the large time behaviour of the interface is given by jX j q ; jxj dq ; i.e., a circular (spherical) steady state. The results of this section provide checks and balances on, and are borne out by, numerical simulations that follow, cf. Figures 1 and 4 and, Figures 7 and 13 for the nutrient-rich and nutrientpoor cases respectively. In particular, p is uniform in space and independent of t when D 0 in the case of unconstrained growth; when part of the tissue boundary impinges on an impermeable surface, say, (representative of a cancer growing against bone for example) then this is no longer the case and nutrient-rich growth with D 0 is no longer shape-preserving, as numerical simulations (not included) confirm.

Weak formulation of the Stokes-flow free boundary problem
Our numerical method for the approximation of the Stokes free boundary problem is based on the following weak formulation of (2.14) and (2.15). We define the function spaces for the L 2 inner product, with the standard modification for vector valued functions.

Numerical scheme
We present a robust numerical method for the approximation of (3.1) based on a finite element discretisation in space. We partition the time interval OE0; T as 0 D t 0 < t 1 < < t N D T , with uniform (for simplicity) timestep D t nC1 t n , n D 0; 1; : : : ; N 1. For n D 0; 1; : : : ; N , we let n h be an approximation to the surface˝.t n / and we set n h D @˝n h , i.e., n h is the boundary of the computational domain˝n h with n h and approximation to .t n /. For each n, we define T n h to be a triangulation of˝n h consisting of the union of d dimensional non-degenerate, non-overlapping closed simplices (triangles for d D 2 and tetrahedra for d D 3) whose intersection is either a shared vertex, a shared complete face (edge) or empty.
We employ Taylor-Hood .P 2 ; P 1 / elements for the Stokes equations and piecewise linear surface finite elements for the approximation of the surface evolution law; to this end we define bulk and surface finite element spaces as follows U n h WD f˚h 2 C.˝n h / W˚hj k 2 P 2 .k/ for each k 2 T n h g d ; P n h WD˚˚h 2 C.˝n h / W˚hj k 2 P 1 .k/ for each k 2 T n h « ; The method we employ is based on parametrising the surface nC1 h over the surface at the previous timestep n h . We will make use of an approximate surface normal n h 2 V n h whose values at the vertices are defined by taking the weighted (by element area) average of the normals to each element in the patch containing the vertex and normalising.
Our numerical method for the approximation of (3.1) for n D 0; 1 : : : ; N reads as follows, given . We update the interior nodes of the triangulation of h by computing the harmonic extension of the boundary displacement (we also tried simply using the computed Stokes velocity to update the interior vertices with no qualitative differences).
As the velocity U h is determined only up to rigid-body motions, we used a null space iterative solution procedure and orthogonalisation of the right-hand side vector with respect to the rigidbody motions, i.e., the three-dimensional in 2d and six-dimensional in 3d spaces of translations and rotations.

Numerical results
We focus here on simple, but representative, initial domains (specifically, ellipses and ellipsoids) in order to exemplify potential modes of instability associated with departures from radial symmetry. A more comprehensive survey of the dependence on the initial data (and encompassing the effects of a tumour abutting a rigid boundary, say, to mimic the presence of bone) would be worthwhile but beyond the scope of the current investigation.

Nutrient-rich case
4.1.1 2d simulations. The free boundary problem in the nutrient-rich case corresponds to (2.14) with 2 R C and q D 0. In all the simulations in 2d we considered an initial domain of the form i.e., a rectangular domain with elliptical caps.
To test the method we first consider the case of zero surface tension, i.e., D 0. As discussed in 2.3, the solutions in this case should correspond to shape-preserving growth. Figure 1 shows the results of a simulation with D 0:5, D 0 and D 1 with the discretisation parameters as given in Table 1. In Figure 1 and all subsequent figures depicting simulation results, we have shaded the interior of the domain according to either the magnitude of the velocity vector or the value of the pressure and shaded the boundary by the magnitude of the mean curvature vector.
We observe, as expected (cf. 2.3.2), shape preserving growth. In order to test the code in the case that the surface tension is nonzero we simulated (2.14) with D 0:5, D 0 and D 2, the domain converging to a disc with area approximately equal to the initial error (< 1% error, results omitted for brevity). The results of a simulation with D 0:5, D 2 and D 1 and the  discretisation parameter values as in Table 1 are shown in Figure 2. We observe marked differences to the shape preserving case, with contraction of the tissue in the curved regions around the elliptical caps and the shape at the end time is close to a disc. Indeed a linear-stability analysis implies convergence towards radial symmetry for > 0; contrasting significantly with the behaviour for D 0. Interestingly, radial symmetry is not in fact attained, for reasons outlined in 5 below. This is evident in Figure 3 where we plot the aspect ratio (AR) of the tissue defined by (using Cartesian coordinates x D .x; y/ T ), against time and observe that the aspect ratio does not converge to 1. Intuition into the reasons for such behaviour can be gleaned by noting that tissue expansion reduces interface curvature and hence the effectiveness of the surface tension term (which in the current context is typically interpreted in terms of cell-cell adhesion).  Table 1. Reading from left to right, snapshots of the solution are shown at t D 0; 0:1; 0:3 and 0:6. We observe growth with the evolution appearing to approach a growing disc.
FIG. 3. Aspect ratio, as defined in (4.2). against time for the Numerical solution of (2.14) in the nutrient-rich case in 2d with D 0:5, D 2; D 1; q D 0 and the discretisation parameter values given in Table 1. We note the aspect ratio converges to approximately 0.95 and that radial symmetry (i.e., an aspect ratio of 1) is not attained, as discussed in 5.   Table 2. The top row shows the surface shaded by the magnitude of H , the middle row a cross section shaded by the magnitude of u and the bottom row a cross-section shaded by p. Within each row, reading from left to right, snapshots of the solution are shown at t D 0; 0:5 and 1. We observe convergence towards a growing sphere with rapid initial contraction near the regions of largest curvature.
simulations we take the initial domain to be an ellipsoid defined by The results are analogous to the 2d simulations with shape preserving growth in the zero surface tension case, Figure 4, and convergence towards (but not attaining) a growing sphere in the positive surface tension case, Figure 5. Figure 6 shows outlines of the surface mesh at the final time of the simulation illustrating the effect of the surface tension term on the shape. The free boundary problem in the nutrient-poor case corresponds to (2.14) with < 0 and q 2 R C . In the first set of simulations in this section we considered the initial domain defined in (4.1), as in 4.1.1. Figure 7 shows the results of a simulation D 0:5, D 0; D 0:5 and q D 1 with the discretisation parameters as given in Table 1. We no longer observe shape preserving growth with the evolution appearing to converge towards a growing disc in agreement with the argument of 2.3.2. The results of a simulation with D 1; D 0:5 and q D 1 and the remaining parameter values as in Table 1 are shown in Figure 8. We observe similar evolution to the case with no surface tension however, in accordance with Section 4.1, we observe contraction of the tissue in the curved regions around the elliptical caps and the convergence towards a growing disc shape is faster than in the case with no surface tension with the domain at the final time significantly smaller in area than  Table 1. Reading from left to right, snapshots of the solution are shown at t D 0; 0:1; 0:3 and 0:6. As in the nutrient-rich case with zero surface tension, the pressure is constant in space and time and takes the value 0:25. We do not observe shape preserving growth, with the growth approaching a growing disc in keeping with the analysis of 2.3.2.
in the case of no surface tension: this is in keeping with the overall area balance in 2d with < 0, where A is the tissue area and L the length of the interface -surface tension acts to reduce L and hence the growth rate of A. As noted in 2.1, the case < 0 is also of interest in this particular limit. In Figure 9 we report on results with D 0:5; D 1; D 0:5 and q D 1, unlike the previous two cases considered in this section we no longer observe growth towards a disc like shape; rather, growth is fastest in the elliptical caps with relatively slow growth in the rectangular region connecting the elliptical caps.
To illustrate the long time convergence towards a stationary disc in the > 0 case and in order to illustrate the instabilities in the < 0 case, we show outlines of the boundary curve shaded by the magnitude of the mean curvature at larger times, namely t D 0; 2; 4; 10 and 20 in Figure 10. We note that whilst numerical instabilities were not observed for any of the other simulations we report on in the manuscript, for this simulation with < 0 we observe the onset of instabilities if we continue the  Table 1. Reading from left to right, snapshots of the solution are shown at t D 0; 0:1; 0:3 and 0:6. We observe growth approaching a growing disc with initial contraction in the elliptical caps. simulation over long times which result in degeneration of the grid and oscillations in the boundary curve such that we cannot continue the simulation (for this example the algorithm broke down at t D 4:4 and this could not be remedied by simply reducing the discretisation parameters, reflecting the presence of instabilities in the model itself as well as in its discretisation). The role played by a negative surface tension coefficient in the stability of the solution and possible regularisations of the model is something we wish to investigate in future work: here we simply emphasise the role of surface compression leading to a buckling-associated fingering mechanism. Figure 11 shows the area of the tissue along with that obtained by solving (4.4) such that where the computed interfacial length was used and the integral was approximated using the Trapezoidal rule. We observe good agreement between the computed area and the 'analytical solution' given by the area law in all three cases.  Table 1. The surface is shown at t D 0; 2; 4; 10; 30 (note that the algorithm failed to converge at t D 4:4 for D 1 hence only three outlines are shown). For > 0 we observe convergence towards a stationary disc, whilst in the case of negative surface tension the onset of a noteworthy buckling-type instability is observed and the final state shown is not a stationary state.
FIG. 11. Computed area (solid lines) and 'analytical area' (red markers) obtained by solving (4.4) with given initial volume and computed interface length for D 0 (cyan curve and red crosses), D 1 (blue curve and red circles) and D 1 (green curve and red diamonds). Good agreement between the solution of the area law and the computed areas is observed in all cases.
As a final numerical simulation in 2d , we consider a circular initial domain corresponding to a disc of radius 1. We considered the cases D 1; 0; 1 with D 0:5, D 0:5 and q D 1 as in the previous example. Figure 12 shows the long time behaviour of the solutions for the three different values of , as expected, radial symmetry is preserved and, in accordance with the area law (4.4), the domains all grow at the same rate independent of the value of . As with the previous example, for < 0 instabilities eventually arise and the algorithm fails to converge for sufficiently large t . This also illustrates that the instabilities seen in Figure 10 are not due to the nonuniform mesh as in this example radial symmetry is preserved until being broken through the onset of the instabilities, the mesh remaining near-uniform until the time the algorithm (and perhaps the model itself) breaks down.  Table 1 for a circular initial domain. The top, middle and bottom rows correspond to D 1; 0; 1 respectively. Snapshots are shown at t D 0; 1; 3; 7; 10. For the D 1 case the algorithm failed to converge at t D 7:3 so only four snapshots are shown: interestingly, near-circular growth is more prone to such breakdown than growth with a preferred direction whereby surface compression can be ameliorated by extension in that direction, as illustrated in Figures 9 and 10. We note that all three cases exhibit similar growth rates in accordance with the area law (4.4). Figures 13-15 show results of simulations in 3d of (2.14) in the nutrientpoor case with D 0:5; D 0:5 and q D 1 for the same cases considered in 2d of zero, positive and negative surface tension. As in the nutrient-rich case, in each of the simulations we take the initial domain to be an ellipsoid defined by (4.3). The numerical parameter values used are given  Table 2. Reading from left to right, snapshots of the solution are shown at t D 0; 0:5 and 1. We observe convergence towards a growing sphere, with rapid initial contraction near the regions of largest curvature.

3d simulations.
in Table 2. The results are similar to the 2d simulations, with growth towards a growing sphere in the zero surface tension case, Figure 13, rapid convergence towards a growing sphere in the positive surface tension case, Figure 14, and elongation in the negative surface tension case with growth fastest in the direction of the longest axis of the ellipsoid, Figure 15. Figure 16 shows outlines of the surface mesh at the final time of the simulation illustrating the effect of the surface tension term  Table 2. Reading from left to right, snapshots of the solution are shown at t D 0; 0:5 and 1. We no longer observe convergence towards a sphere, with the aspect ratio of the final shape larger than that of the initial ellipsoid. on the shape. In 3d dV dt D qA C V; (4.5) with < 0, where V is the tissue volume and A the interfacial surface area: in this case we highlight that surface compression . < 0/ drives an increase in A and hence an increase in the overall growth rate.

Discussion
Tissue growth motivates the investigation of moving boundary problems that generalise classical formulations. We have explored two such generalisations of established Stokes-flow frameworks.
The nutrient-rich model contains a single additional term, namely a uniform volumetric source, so that dA dt D A; dV dt D V; (5.1) in 2d and 3d , where A denotes the area of the tissue (in 2d ) and V the volume of the tissue (in 3d ), with benign modifications to the classical case beyond the exponential growth exhibited in (5.1). Indeed, defining X ; U ; P and˙by x D e as in 2.3.2, leads in the nutrient-rich case to the classical Stokes free boundary problem (i.e., that without a volumetric source, cf. the analysis of the mathematically related problem in [5]), so that the two-dimensional case is amenable to complex-variable methods (see [6] and references therein). It is important to note that remains finite as t ! 1 so that so that the large behaviour (i.e., radial symmetry) is not fully attained. In the nutrient-poor case < 0 the same change of variables eliminates the volumetric sink but preserves the constant surface source. In this case ! 1 as t ! 1; the presence of the surface source would seem to militate against the application of complex-variable methods in the two-dimensional case, though further consideration of such approaches would be in order.
Of the two cases considered here, the nutrient-poor one is undoubtedly the more challenging, due both to the competition between the surface source and volumetric sink made explicit by (4.5) and to the relevance of surface compression. The latter drives the flow away from radial symmetry, with both the numerical results and a preliminary analysis suggesting the scope for parallel-sided growth away from the tips. In any case, the < 0 case is expected to be mathematically ill-posed but, informally, not badly so; in consequence it raises a number of worthwhile questions both in its own right and in terms of the identification and analysis of suitable regularisations.
Instabilities associated with tumour surface morphologies are often viewed as indicative of enhanced invasiveness and hence of a worse prognosis. Such effects warrant further investigation, though caution should naturally be exercised in drawing detailed conclusions from a mathematical framework as simple as that described above (to cite just one example, account should also be taken of two-phase effects, whereby the growing tumour is surrounded by healthy tissue). Nevertheless, the specific buckling-associated instability mechanism evidenced above should be of fairly general relevance and the results should be of most immediate specific applicability to the in-vitro growth of tumour spheroids and embryoid bodies, for example (albeit with similar provisos regarding the maximal simplicity we have sought in the model frameworks treated above).
The numerical simulations we present are consistent with the analytical results of 2.3 as well as the area and volume laws stated in (5.1). Further investigation of numerical methods adapted to evolution laws with negative surface tension is warranted as this is much less studied than the positive surface tension case. The approximation of higher-order differential operators on the surface is also of interest as it is expected to arise when we consider regularised models for the case of negative surface tension.