The buckling of capillaries in solid tumours

We develop a model of the buckling (both planar and axial) of capillaries in cancer tumours, using nonlinear solid mechanics. The compressive stress in the tumour interstitium is modelled as a consequence of the rapid proliferation of the tumour cells, using a multiplicative decomposition of the deformation gradient. In turn, the tumour cell proliferation is determined by the oxygen concentration (which is governed by the diffusion equation) and the solid stress. We apply a linear stability analysis to determine the onset of mechanical instability, and the Liapunov–Schmidt reduction to determine the postbuckling behaviour. We find that planar modes usually go unstable before axial modes, so that our model can explain the buckling of capillaries, but not as easily their tortuosity. We also find that the inclusion of anisotropic growth in our model can substantially affect the onset of buckling. Anisotropic growth also results in a feedback effect that substantially affects the magnitude of the buckle.


Introduction
It has been hypothesized, for a long time, that intratumoural blood vessel abnormalities can impair blood flow into the tumour, and consequently decrease the effectiveness of chemotherapeutic treatment (Jain 2008). Tumour blood vessels can have many types of abnormality, including chaotic direction and entanglement, a high level of branching, overly porous capillary walls and chaotic blood flow. One particular type of abnormality is vascular buckling and collapse. In this paper, we develop a mathematical model of this phenomenon. Vascular buckling in the circumferential direction has been well-documented by experimentalists. Goldacre & Sylven (1962) and Boucher & Jain (1992) conjecture that such planar buckling is caused by compressive stress arising from tumour cell growth. Indeed Griffon-Etienne et al. (1999) and Padera et al. (2004) have observed that drug-induced apoptosis around a collapsed lumen caused the lumen to reopen. Vascular buckling and collapse has been modelled in two major ways. In the first way, the capillary wall is modelled as a thin shell subject to an external hydrostatic pressure (Stoll 2003). In the second way, Araujo & McElwain (2004) model the stresses accumulating in a growing tumour and conjecture that buckling occurs when the stress exceeds a certain threshold.
The high degree of tortuosity of the tumour vasculature has also been well-documented by experimentalists. Nagy et al. (2009) conjecture that this is primarily due to excessive lengthening of the growing capillary, so that if the ends of the capillary are fixed, it is forced to coil to accommodate the extra length. In this paper, we build a model that seeks to explain vascular buckling (both planar and axial) as a mechanical instability resulting from the rapidly proliferating tumour cells. We wish to investigate whether differential growth induced by the oxygen profile can explain capillary buckling. It has already been shown that differential growth can explain the buckling of spherical shells in Ben Amar & Goriely (2005), and the wrinkling of tubular mucous-lined cavities in Moulton & Goriely (2011) and Li et al. (2011). We neglect the capillary wall as a distinct structure, although we include it in a subsequent paper. This is because we wish to investigate the hypothesis of Fung (2004) that the stability of blood vessels derives more from the surrounding tissue than the strength of the vessel wall. Our model will seek only to explain vascular buckling in terms of tumour cell proliferation (and the resulting strain and stress) in the neighbourhood of a single capillary.

Model
We model the buckling of capillaries embedded within a larger tumour mass with many other capillaries. However, we include only a single capillary, and its immediate neighbourhood, in our model. (Some researchers, such as Bertuzzi et al. (2003), would call this a tumour cord, but we do not use this term in order that the language is as transparent as possible.) Thus we consider the tumour to be a cylindrical annulus, with inner radius R 1 and outer radius R 2 . In the axial direction, the boundary of the annulus is given by Z = ± 1 2 Z 0 . Thus, we consider R 2 to be approximately half the distance of our capillary to its nearest neighbour. In the results, particularly §4, we will discuss the conditions under which our restriction of the geometry to that outlined above is accurate. Because we are neglecting the capillary wall as a distinct mechanical structure, the 'capillary' is merely an annular hole in the material. We denote the reference configuration of the tumour body by V 0 , and the current configuration by V. The geometry of our model is shown in figure 1.

(a) The mechanical structure of tumours
There are two major types of continuum model of the mechanical structure of tumours: multi-phase models and single-phase elastic models. Multi-phase models decompose the tumour into separate phases: typically the extracellular matrix, cells and interstitial fluid (Roose et al. 2003;Chaplain et al. 2006), whereas single-phase elastic models approximate the tumour as a homogenous nonlinear Here, R ∈ [R 1 , R 2 ] is the tumour annulus, where R 1 is the radius of the capillary wall, and R 2 is the outer radius of the annulus (of the order of half the mean intercapillary distance). We neglect the capillary wall as a distinct capillary structure. Z 0 is of the order of the mean interbranching distance of capillaries. elastic material. We choose the latter for our model, chiefly because most of the major strengths of the multi-phase model are not realized when the tumour is in quasi-static equilibrium over the growth timescale (which is a matter of days). In particular, we neglect poroelastic effects from the model because the hydrostatic pressure gradient and fluid flux in tumours is often not nearly as great as in normal tissue, mostly because the lymphatic system of tumours often functions poorly (Boucher & Jain 1992). In addition, the poroelastic relaxation timescale is a matter of seconds (Netti et al. 2000), so that any local fluid pressure gradients resulting from the growth quickly dissipate over the growth timescale. The viscoelastic relaxation timescale is similarly short.
Let u be the displacement, and F = I + V 0 u be the deformation gradient. Note that, in this paper, the components of the vectors and tensors at a point with reference coordinates (R, Q, Z ) are always relative to the cylindrical orthonormal basis vectors {e R , e Q , e Z }. We model the tumour as a Blatz-Ko (BK) material, with strain energy function parameters were fitted using data from confined compressions tests on tumour spheroids by Helmlinger et al. (1997), Netti et al. (2000) and Roose et al. (2003). These experiments revealed that the single-phase elastic model is most accurate for collagenous tumours where the extracellular matrix does most of the loadbearing, particularly glioblastomas. It is important to note that the BK model is compressible. The confined compressions experiments demonstrated that some tumours are compressible owing to the drainage of the fluid from the pores. It is these experiments that motivated our choice of the BK material, as it was originally developed to model porous, rubber-like materials.

(b) Growth distribution
There are many variables that influence the tumour growth. These include nutrients, such as oxygen and glucose, chemical growth factors and solid stress (Helmlinger et al. 1997). In our model, we will simplify the rate of growth to be dependent on the concentration of only one nutrient (oxygen). There is considerable variation in the oxygenation of tumours, probably due to the heterogeneity of the vasculature (Vaupel et al. 1989).
Let r g be the mass growth, per unit of current mass, in the current configuration. We will determine constitutive relations for r g in terms of the concentration of the oxygen, and then use r g to determine the growth tensor. However, before doing this, we must first model the distribution of oxygen.

(i) Mass growth and nutrient diffusion
We denote the concentration of oxygen per unit current volume by c, and model it using a reaction-diffusion equation over current coordinates where G n is the rate of uptake of oxygen per volume of tumour, and D is the diffusion coefficient. Here, and in subsequent expressions, V is the nabla over the spatial coordinates (and V 0 is the nabla over the reference coordinates). We have already assumed in §2a that fluid flux is negligible. Thus, the distribution of oxygen is diffusion-dominated, rather than advection-dominated (as in ordinary tissue). We may model the nutrient distribution to be in quasi-steady equilibrium, because the diffusive timescale is much less than the growth timescale, which we justify as follows. The diffusive timescale [t] d = [R] 2 /D gives the timescale over which diffusion is significant. The diffusion coefficient is approximately 10 −9 m 2 s −1 (Roose et al. 2003). The lengthscale is approximately [R] = 10 −4 m (Konerding et al. 2001). Hence, we find the diffusive timescale over this distance is of the order of 10 s. Because the growth timescale is a matter of days (Casciari et al. 1992), we may model the oxygen concentration as having reached a quasi-steady-state equilibrium, i.e.
DV 2 c − G n = 0 over V. (2.3) We note that the nutrient concentration profile varies with respect to time as the tumour grows owing to the deformation of the geometry. However, in light of the above analysis, we assume that the nutrient levels quickly equilibrate. We set c at the inner boundary (with reference coordinate R = R 1 ) to be the blood oxygen concentration, and the flux of c at the outer boundary (with reference coordinate R = R 2 ) to be zero, which translates to the following boundary conditions for c : where n is the unit outer normal in spatial coordinates. We model the rate of uptake of oxygen per volume of interstitium using Michaelis-Menten kinetics. Justification for this model can be found in the experiments of Casciari et al. (1992), who demonstrated that the uptake per volume depends approximately linearly on the oxygen concentration for hypoxic spheroids, but that it asymptotes to a limit as c → ∞. We thus find G n = l n c(k n + c) −1 , where l n and k n are constants. The rate of growth per mass is also modelled using Michaelis-Menten type kinetics, using parameters from the same experiments. Hence, our expression for the rate of growth is r g = (l c c)(k c + c) −1 , where l c and k c are constants. Note that the experiments of Casciari et al. (1992) give the net rate of growth of the tumour cells, which means we do not need to include any terms accounting for cell death, because l c accounts for it already. It needs to be noted that we have neglected the effect of solid stress on the total mass growth. The evidence for this is mostly phenomenological, and it is very difficult to construct a mathematical model from the data (Helmlinger et al. 1997). Including this would unnecessarily complicate our model. However, we will, further below, consider the effect of stress on the direction of mass growth.

(ii) Growth tensor
In this section, we derive constitutive relations for the stress in terms of the mass growth r g . The deformation gradient is multiplicatively decomposed into an elastic deformation F e and a growth tensor F g , i.e. F = F e · F g . (2.5) The rationale for this ansatz is that the growth tensor causes the material to pointwise expand, resulting in a hypothetical unstressed intermediate state V h . The elastic deformation then corrects the growth tensor to ensure the integrity of the material, but internal stress buildup results as a consequence of this correction. In other words, we substitute F e into the elastic constitutive relations rather than F. This approach was first used to model volumetric growth by Rodriguez et al. (1994), and has recently been applied to tumoural growth by Ambrosi and co-workers (Ambrosi & Mollica 2001;Ambrosi & Guana 2005). It has been used to model the buckling of growing cylindrical objects such as ours by (among others) Li et al. (2011) and Moulton & Goriely (2011). The growth from V 0 to V h is assumed to be volumetric, such that the density of the grown material is the same as the density of the material in the reference configuration. Thus, if j g and j are, respectively, the strain energy per unit volume of V h and the strain energy per unit volume of V 0 , then j = hj g , where h = det(F g ). To see why this is the case, one can imagine fixing a small segment of material dm in the reference coordinates. The volume of dm grows by a factor of h from V → V 0 , but the strain energy in dm remains the same. It follows from the above constitutive assumptions that j g = j BK (F T e · F e ), where j BK is the strain energy function of the BK material given in (2.1). We obtain the first Piola-Kirchhoff stress tensor P by differentiating the strain energy per reference volume j with respect to the deformation gradient, i.e. P = h(vj g /vF). We find

(iii) Growth tensor and mass growth
It remains for us to determine a constitutive relation between F g and r g . Throughout this section, we determine constitutive relations for the variables at a point with fixed reference coordinate X . It follows from the definition of r g that d dt (rJ ) = rJr g , where r is the density. (2.7) We assumed earlier that the growth tensor from V to V h is purely volumetric, so that the rate of change of the density over V h is zero, i.e. d/dt(rJ e ) = 0, where J e = det(F e ). Because J = hJ e , (2.7) may be simplified to Prior to buckling, the deformation will be assumed to be purely radial. This means that, prior to buckling, the stress tensor is diagonal. We assume that the growth tensor is also diagonal, with non-zero elements z R = F g,RR , z Q = F g,QQ and z Z = F g,ZZ . Hence, h = z R z Q z Z and (2.8) becomes This motivates us to define the anisotropic growth multipliers for which we outline constitutive relations in the next section. On integration of (2.10), we find We have assumed z a = 1 when t = 0.
(iv) Anisotropic growth In this section, we determine constitutive relations for the anisotropic growth multipliers. These govern the responsiveness of the direction of growth to the prevailing stress field. There is abundant experimental evidence that growth occurs preferentially in the direction of least stress (Helmlinger et al. 1997;Chen et al. 2004). However, there is a paucity of quantified experimental data from which we may construct a model. We therefore use simple constitutive relations of the form where a ∈ {R, Q, Z }, T is the Cauchy stress tensor (which is diagonal prior to buckling), K t is the bulk modulus of the tumour and A is an undetermined parameter. The parameter A represents the responsiveness of growth to the prevailing stress field. If A = 0, the growth is isotropic, whereas in the limit as A → ∞ all of the growth occurs in the direction of least compression. These constitutive relations are continuously differentiable and satisfy G R + G Q + G Z = 1, because, as stated previously, we assume that the stress does not affect the net level of growth.

(c) Mechanical equilibrium and boundary conditions
The system is in mechanical equilibrium, so that over V 0 (2.14) We stipulate the axial displacement to be zero at Z = ± 1 2 Z 0 , i.e. u Z = 0. (2.15) We have already stated that the top and bottom faces of the annulus correspond to branch points of the capillary, so that (2.15) is equivalent to stipulating that the axial displacement at the branch points is zero. We ignore any non-radially symmetric or torsional effects of the branch on the growth/stress in our model. Because the interstitial fluid pressure is approximately equal to the microvascular fluid pressure (Boucher & Jain 1992), we assume that the mechanical stress at the inner boundary of the capillary wall is zero, i.e. P · N = 0 at R = R 1 , where N is the unit outer normal. (2.16) The boundary condition at R 2 is null-displacement, i.e. (2.17)

Buckling analysis
We will investigate the mechanical stability of our system according to the energy criterion of Lagrange. According to this energy criterion, a system is stable if the energy functional has a local minimum there (Marsden & Hughes 1983). We assume that initially the deformation and nutrient distribution are purely radial, and that the growth tensor is diagonal with respect to the cylindrical orthonormal basis. We denote this state I. The primary (growth-induced) stresses and strains of state I are in the plane perpendicular to the cylinder's axis, so that our model differs from the classical problem of the Euler buckling of a column. The standard cylindrical orthonormal basis of the reference configuration may also be naturally used as a basis for state I because the displacement is purely radial. Thus, all of the following vectors and tensors are with respect to this basis. We investigate the mechanical stability of state I by performing a perturbation expansion for a small increment in the displacement about state I, in the manner of Ogden (1984). Upon doing this, we obtain a second state, which we denote as state II, i.e.
Here, u 0 is the displacement of state I, and W 0 is the total mechanical energy of state I. The first-order increment in the displacement (u 1 ) must be admissible: meaning that it is in H 1 (V 0 ) and satisfies the following boundary conditions: and We note that here, and after, V 0 u 1kl is the (k, l) component of the second-order tensor V 0 u 1 . The incremental stress tensors vP/vF are, in the terminology of Ogden (1984), the fixed-reference elastic moduli. The derivative is performed in state I. In our investigation of the onset of buckling, we do not consider a small increment in the nutrient distribution c or mass growth r g . This is because the elastic deformation timescale (which is assumed to be much less than a second) is much more rapid than the diffusive timescale (which we have already seen is a matter of seconds) and the growth timescale (a matter of days).
Because state I is in mechanical equilibrium, W 1 is identically zero. The second term W 2 is quadratic in u 1 . We use the Trefftz criterion to determine the onset of buckling (this is discussed in detail in Ogden (1984), although it is not explicitly identified). This states that buckling occurs when, for all admissible v, the first variation of The Euler-Lagrange equation that results from (3.6) is In other words, the existence of a solution u 1 to (3.7) and (3.3)-(3.5) indicates the destabilization of state I. We note that, if the solution u 1 exists, then W 2 (u 1 ) = 0 (which can be seen from substituting v = u 1 into (3.6)).

(a) Fourier decomposition
We investigate the stability of the system to buckles of particular wavelengths by assuming u 1 = au (n,u) for some a ∈ R, where and u (n,u) = 1 at R = R 1 . We are here defining U n,u to be the set of all displacements that can be written in the form (3.8). We justify the above assumption on the structure of u 1 in the electronic supplementary material, §S1. When such a solution u (n,u) exists, we denote it an (n, u) eigenmode. We find the eigenmode u (n,u) by gradually incrementing the critical parameter (the time) until there exists a solution of the requisite form to (3.7) and (3.3)-(3.5). We denote the critical time at which such an eigenmode exists by t (n,u) . We are interested in determining the (n, u) eigenmode even if this is not the first to go unstable. There are two fundamental reasons we do this. First, it will be seen in our results that often the modes of very fine wavelength go unstable first. However, we do not expect our continuum approximation to be accurate for such fine wavelengths. In particular, we expect the cellular phase to dispel such fine wavelengths (as the cellular radius is approximately one-third to a half the capillary radius). Second, there may be small inhomogeneities which bias the system towards some modes more than other modes.
Note that if (n, u) is not the first mode to go unstable, then u (n,u) does not minimize W 2 , but is a saddle point. However, u (n,u) does minimize W 2 with respect to the set of all displacements in U n,u . It can be shown that u (n,u) is infinitely continuous and differentiable. This essentially follows from the fact that (3.7) and (3.3)-(3.5) transforms into an elliptic ordinary differential equation in R under the Fourier decomposition in (3.8). The ellipticity follows from the fact that, as long as q < −1 (which holds for all the tumours we apply our model to, as can be seen from the data in the electronic supplementary material, §S1), there exists a positive K such that for a ∈ {R, Q, Z }, This ellipticity condition (3.9) may be inferred from the expression for the incremental stress tensor outlined in the electronic supplementary material, §S1. Our method of numerical solution is outlined in the electronic supplementary material, §S7.

(b) Magnitude and stability of the buckle
We determine the leading order of the magnitude and stability of the buckle by using a weakly nonlinear analysis about the critical time of buckling. Owing to a shortage of space, we provide only a summary of our method here: we outline our method more fully in the electronic supplementary material, § §S4-S6. In essence, we require the displacement to minimize the higher-order terms of the energy, a method that has been used in the context of shell theory by Koiter (2009).
We stipulate the leading order of the increment in the time to be O(e 2 ), i.e. t = t (n,u) + 1 2 e 2 t 2 + · · · , which is necessary for the perturbation expansion to be well-matched. We find that W 1 is identically zero as a consequence of state I being in mechanical equilibrium, W 2 is identically zero as a consequence of the Trefftz condition and W 3 is identically zero as a result of the regularity of the geometry. We determine the magnitude of the buckle (i.e. a) from W 4 . We do this by using the Liapunov-Schmidt reduction. We first fix a and require that u 2 renders W 4 stationary. The resulting Euler-Lagrange equation allows us to determine u 2 as a function of a and t 2 . We may then substitute this expression for u 2 back into W 4 , obtaining an expression (which is independent of u 2 ) of the form where {C 1 , C 2 , C z 3 , C n 3 } are constants given in the electronic supplementary material, §S4. The constants C z 3 and C n 3 are due to the feedback of the buckle on the growth tensor, which we discuss further in §5c. We stipulate that the system enters the buckled state at kt 2 , for some constant 0 ≤ k ≤ 1. The constant k is undetermined by our weakly nonlinear analysis. We determine a by requiring the derivative of W 4 with respect to a to be zero. On taking this derivative, the integral t 2 kt 2 ad(t 2 ) is taken to be independent of a, because this represents the cumulative effect of the displacement on the strain energy in the past, but we only require the displacement to minimize the strain energy in the present. We obtain two non-zero solutions of the form (3.11) The buckled equilibrium states corresponding to these two solutions are identical: one may be obtained from the other by rotating by p/n. The stability of the buckled equilibrium state is governed by the second derivative of W 4 with respect to a, i.e. v 2 W 4 va 2 = 3a 2 C 1 + C 2 t 2 . (3.12)

Purely radial solution
In §4, we investigate the solution for state I predicted by our model (i.e. the purely radial initial state, which we indicate with the index 0). We must solve state I numerically, using the method outlined in the electronic supplementary material. We fix the parameters at typical values, as outlined in table 1. Because we are modelling an in vivo phenomenon, there is a paucity of data against which we may evaluate our models. However, the following criteria are of use in guiding our evaluation. It is important to note that we cannot prove any of these criteria, Table 1. The default parameter values we used in our results, unless otherwise indicated. We use R 1 as our lengthscale and l −1 c as our timescale. We do not need to explicitly state a value of l −1 c because the timescale factors out, although we note that l −1 c varies from 10 h to 600 h (Freyer & Sutherland 1985). The shear modulus G t similarly factors out. Refer to the electronic supplementary material, §S8 for an explanation of these parameters and the papers from which they are sourced.
parameter value description R 2 10 −4 m outer radius of annulus (for finite geometries) R 1 1.5 × 10 −5 m inner radius of annulus q −7.9 tumour elastic parameter f 0.2 tumour elastic parameter E t 7800 Pa Young's modulus of tumour l n 0.01 mMs −1 rate of uptake of oxygen k n 4.6 × 10 −3 mM Michaelis constant for oxygen uptake k c 8.3 × 10 −3 mM Michaelis constant for cell growth D 1.5 × 10 −9 m 2 s −1 diffusion coefficient of oxygen in tumour c i 0.1 mM oxygen concentration at R 1 and indeed occasionally some of them may not hold. However, we expect them to hold in most circumstances.
(1) It must be possible to accurately re-integrate the predictions of our (single capillary) model into a large-scale tumour model. (2) We expect the displacement u 0 to be negative at R 1 .
(3) The displacement u 0 should be positive for large R.
We justify these criteria in the electronic supplementary material, §S3. In particular, we explain there what we mean by 're-integration' and 'large-scale tumour model' in criterion (1).

(a) Results
We plot the solutions (for the purely radial state I) of the single-capillary model, as explained at the beginning of this section, for a typical selection of parameters in figure 2. These results reveal numerous problems; most importantly, they reveal that the results are not well-integrated (i.e. criterion (1) is violated). We discuss this in more detail. Suppose that we wished to integrate the above results into the annuli surrounding two adjacent capillaries, which we denote by C A and C B . Firstly observe from figure 2 that the stress grows very rapidly as time progresses, chiefly because the outer boundary condition is fixed but the growth is substantial (the growth tensor at R 2 at t = 1 is approx. 1.3 × I ). Because the radial pressure at R 2 is large, one would expect in most circumstances that in the large-scale tumour model, the compressive stress at the intersection of the annuli surrounding C A and C B would push them apart. This contradicts the null-displacement model, which must insist that the intercapillary distance remains fixed, and it is why we expect the intercapillary distance to increase (i.e. criterion (2)) in most circumstances. Another reason why we think criterion (1) is not satisfied is that we expect the stress, strain and growth to still be substantial at points outside the annular geometry: there is no good reason to think that this would not have a substantial effect on the solution within the annulus. More generally, the restriction of the geometry to an annulus is highly vulnerable to non-annular inhomogeneities in the variables. Because the stress is a nonlinear function of both the growth and strain, and because both the strain and growth at the outer boundary R 2 are already substantial, we expect that inhomogeneities in these variables (owing to factors such as some points on the outer boundary of C A being closer to C B than other points) could have a very large impact on the results. Finally, in our buckled results further below, we find that most of the modes go unstable before t = 0.5, which seems too rapid, because this is less than the characteristic time of cellular proliferation (which is our timescale). It was demonstrated in Konerding et al. (2001) that there is huge variation in the intercapillary distances of different tumours (and also within tumours). We are therefore interested in investigating the behaviour of our model as R 2 /R 1 tends to infinity. We denote the model with R 2 = ∞ the 'infinite-geometry' model. This model is of interest because there is no 'confinement' owing to a fixed outer boundary condition, for which we have argued in the previous paragraph there are problems (when R 2 /R 1 is approx. 6.67). Indeed in the infinite-geometry model buckling may occur only owing to the stress resulting from the differential growth profile. Furthermore, it may be seen from figure 3 that the displacement is positive for large R, so that the infinite-geometry model satisfies criterion (3) (unlike the finite-geometry model). We expect that, as long as the distance between capillaries is sufficiently large, the infinite-geometry model is well-integrable (of course, when 'integrating' the results according to the method in the electronic supplementary material, §S3, one would only keep the [R 1 , R i ] annular subset of the infinite-geometry solution). Of course, we do not expect this to be true if the intercapillary distance is small. A weakness of the infinite-geometry model is that the displacement at R 1 is positive, contradicting criterion (2). It is not surprising that our model predicts this as such, because the material needs to expand at   Figure 4. In the above simulations, we varied the growth anisotropy parameter A from 0 to 50 (with increment 1) and plotted data at t = 1 (relative to the characteristic cell proliferation timescale). We used the infinite-geometry model. The growth is isotropic when A = 0 and becomes increasingly responsive to the prevailing stress field as A increases. At the top is the displacement and the displacement derivative at R 1 (the radius of the cavity). At the bottom is the growth tensor and azimuthal stress at R 1 . Observe that as A gets larger, the differences in growth become more pronounced and the displacement at R 1 is more negative. (Online version in colour.) R = R 1 to relieve the enormous hoop stress. In figure 4, we investigate whether the inclusion of growth anisotropy in our model (which we have noted is very well-attested experimentally) can correct this problem. It can be observed that as A increases (and the growth thereby is more anisotropic), the displacement at R 1 switches from positive to negative, as we desired. We have already noted that it is very difficult to parametrize A from the literature; so we estimate A = 10 from these results. We have chosen 10 as the approximate value for which the displacement at R 1 switches from positive to negative at t = 1.

Buckled results
In this section, we outline our buckled results. There are two major emphases of our results. First, we will investigate the geometry of the buckling, i.e. the nature of the buckled profiles, and the relative stability of the different eigenmodes. Second, we investigate the effect of growth anisotropy on buckling, with a particular emphasis on the feedback effect of the buckle on the growth. In a subsequent paper, we will explore the implications of our model for therapy. Typical parameter values are outlined in table 1. The magnitude and stability of the buckled state is determined, using the method summarized in §3b. We expect considerable variance in the axial wavenumber u. This is because u is of the order of 2 mp/Z 0 , where Z 0 is the interbranch distance and m is a positive integer. Using our data in electronic supplementary material, §S8 we find the interbranch distance is of the order of 130 mm, but varies considerably. We take m = 1 to find u is approximately 1 (note that we have non-dimensionalized such that R 1 = 1).
We briefly comment on the asymptotic convergence of the finite-geometry model results to the infinite-geometry model results, to continue the discussion in §4a. The results are plotted in figure 5. It can be seen that the convergence is slow: for example, for the critical times to be within 10 per cent of their asymptotic limit, the outer radius R 2 must be greater than a millimetre, which is approximately 10 times the mean outer radius in table 1. However, as we have already discussed in the previous section, for large R, the infinite-geometry model is superior in some respects to the finite-geometry model. In the rest of this paper, we will study only the infinite-geometry model.

(a) Buckled profiles
We produce figures of the buckled shape for a typical selection of parameters (given in table 1). The shapes of the planar buckles of modes 2-5 are plotted in figure 6. We observe that as the mode number increases, the amplitude of the buckle gets smaller. Indeed, it may be shown using a boundary layer analysis (as outlined in MacLaurin et al. (in preparation)) that the asymptotic order of the amplitude a is n −3/2 √ t 2 as n → ∞ (for small t 2 ). It can also be observed that the decay as R → ∞ of the buckle is much slower for low n than high n, which is also predicted by our boundary layer analysis. Diagrams of the profiles of axial buckles are plotted in figure 7.

(b) The preferred buckling modes of the system
We now discuss the preferred buckling modes of the system. It may first be observed in all our results that the plotted modes go unstable at times of the order of two to 10 times the growth proliferation timescale, which is in good broad agreement with the experimental literature (Padera et al. 2004). In the following discussion, we restrict ourselves to the situation where A is sufficiently great (i.e. 10 or more) that u 0 (R 1 ) is negative at the time of buckling (which we consider to be more realistic, as stated in criterion (2) of §4).
We do not expect the system to necessarily adopt the first mode to go unstable. This is primarily because the wavelengths of the first modes to go unstable may be so fine that subcellular processes dominate, and our model is not  Figure 5. The effect of varying the outer radius R 2 from 50 mm to 5 mm on the buckling of planar wavenumbers n ∈ {2, 3, 4, 5, 10, 20}. The growth is isotropic. We plot the critical time of buckling in the top plot, and the magnitude a of the buckle in the plot on the left (we plot a 2 n 3 /t 2 because we know from our boundary layer analysis that this asymptotes to a finite constant as n → ∞). In these two plots, we plot the asymptotic limit of these variables as the wavenumber n → ∞ in a black 'x'. In the critical time plot, the asymptotic limit as R 2 → ∞ for each of the modes is plotted in a dashed line. In the postbuckled stability plot on the bottom right, we plot an 'x' if stable buckled states exist, and a 'o' otherwise. The dimensional parameter values are given in table 1. (Online version in colour.) physically accurate. The second reason is that we expect there may be physical inhomogeneities that bias the system more towards one mode than another. Furthermore, as will be seen below, sometimes the energy relaxation afforded by modes that buckle later is greater (to leading order in e) than the energy relaxation afforded by the modes which buckle first. We therefore cannot be certain which modes are preferred by the system, and we are still interested in studying a mode even if it is not the first to destabilize.
It can be observed in the results of figures 8 and 10 that the critical times of the modes are often highly clustered. These results contrast the modelling of capillary buckling in which the capillary is a thin ring: in the latter, the critical pressure at which the ring buckles is proportional to the square of the wavenumber (Stoll 2003). In fact, the key reason for this difference is that the thin ring model cannot accommodate the development of a boundary layer. The clustering of the modes in our model was investigated using a boundary layer analysis, using the method outlined in MacLaurin et al. (in preparation). That is, we fix n = b 1 2 and u = b 2 2, and consider the behaviour of the system as 2 → ∞. It can be observed in our results that for almost all values of b 1 and b 2 , the critical times asymptote to a finite limit. In other words, as the time approaches the asymptotic limit, the system becomes unstable to wrinkles of a progressively finer wavelength. Sometimes, this limit is the infimum of all the critical times: for example, this occurs for large A in the planar buckling Figure 6. The planar buckled profiles. We plot the displacement up to quadratic accuracy, i.e. u = u 0 + eau (n,u) + 1 2 e 2 a 2 u 2,a + t 2 1 2 e 2 u 2,p and t = t (n,u) + 1 2 e 2 t 2 + · · · . The eigenmode (i.e. u (n,u) ) of (a) is mode 2; that of (b) is 3; that of (c) is 4 and that of (d) is 5. We implemented the infinitegeometry model, with anisotropic growth and A = 10 (the data derive from the A = 10 results in figure 8). We have chosen e and t 2 such that 1 2 e 2 t 2 = 0.5, and k = 1. These values are perhaps slightly beyond the regime of accuracy of our perturbation expansion: they were chosen so that the buckled profile may be more easily visualized. The lines represent contours of constant radius (in the reference coordinates) spaced evenly from R 1 to 1.05 × R 1 . The parameters are given in table 1. •(1,1) •(1,10) •(1,0.1) (1,0.1) (1,1) (2,1) (3,1) Figure 10. The effect of varying the growth anisotropy parameter A on axial buckling. We plot the critical times of buckles of mode (n, u). Here ∞(b 1 , b 2 ) is the asymptotic limit of the critical time of (2b 1 , 2b 2 ) as 2 → ∞. The parameters are given in table 1. (Online version in colour.) in figure 8. This behaviour is analogous to the development of surface waves in a compressed infinite elastic half-space (Biot 1966). In that case, if there is a homogeneous uniaxial compression in an infinite elastic half-space, then all the wavelengths (for transverse buckles) become unstable simultaneously because there is no natural length scale and therefore the wavelengths are indistinguishable. Dervaux & Ben Amar (2011) similarly observed a vanishing of the wavelength in the growth of a soft ring surrounding a stiff substrate. Of course, we do not expect infinitely fine wrinkling to actually occur in vivo; largely because our approximation of the growing tumour as a continuum is not accurate for very fine wavelengths (i.e. less than the radius of a typical cancer cell). We hypothesize, using our results, that the system might prefer buckles with a small planar mode, for the following reasons. First, it can be observed that, after restricting ourselves to the case where u 0 (R 1 ) is negative, the order of buckling of the planar modes (in figure 8) is not affected by A. Similarly, the order of buckling of most of the axial buckles is also not affected by variance of A, as long as u 0 (R 1 ) is negative. In both the planar and axial buckling, we typically (but not always) find U n,u goes unstable after U m,u if n < m, which might seem to suggest that the system prefers buckles of short planar wavelength to buckles of long planar wavelength. Moulton & Goriely (2011) also observed that buckles of short planar wavelength tended to go unstable last in their bifurcation analysis of an unconstrained (and unpressurized) cylindrical tube growing anisotropically (furthermore, their results seem to strongly suggest the existence of a boundary layer for large n). Similarly, Li et al. (2011), in their model of the wrinkling of a mucosal tube surrounded by a stiff outer layer, found that the system tended to prefer planar buckles of relatively high mode (approx. 6). However, as explained above, we do not expect our approximation of the tumour as a homogeneous material to be accurate for buckles of very small wavelength. Furthermore, the leading order of the decrease in energy induced by the system moving from state I to the buckled state (i.e. W 4 ) is almost always greater for buckles of type (n, u). This means that, when the critical times of the buckles are sufficiently clustered that we require only the leading order of the energy W 4 to accurately compare them, we almost always find that the (n, u) buckles are energetically more favourable. Thus, for all of the reasons above, we tentatively expect that the system prefers buckles with a small planar mode (n). This is consistent with the fact that Stoll (2003) observed buckles only with n = 2 and n = 3.
We cannot conclude with certainty from these results whether growth contributes to capillary tortuosity; however, it seems unlikely, for the following reasons. The modes that are most likely to contribute to capillary tortuosity are the ones with n = 1, because this is the only planar wavenumber for which the centre of each planar buckled slice is not on the Z axis. However, as we argued above, modes of this form (i.e. of the form (1, u)) are usually the very last to go unstable. Indeed, over the time range of our search (i.e. with t < 10), if the system goes unstable to a buckle of this form, then the displacement u 0 at R 1 at the time of buckling is positive, contradicting criterion (2). One might conclude, using the argument of the previous paragraph, that these buckles might be energetically more favourable. However, this argument is valid only for buckles for which the critical times are closely clustered together (so that one requires only W 4 to accurately compare the energies of the buckled states). However, buckles of this form are almost never closely clustered with buckles for which n = 1. This can be seen from the fact that, even when buckles of this form are observed (for low A), the growth h at R 1 at the critical time of buckling is typically at least 100 times as great as the growth at R 1 at the critical time of the planar buckles. In conclusion, it seems unlikely that tumour growth contributes significantly to capillary tortuosity, particularly because there are other plausible explanations for this phenomenon. One possible explanation is that in Nagy et al. (2009), i.e. that the capillary coils due to excessive growth in the axial direction.

(c) Growth anisotropy and feedback effects on buckling
The response of the axial buckling to variation of the growth anisotropy parameter A is plotted in figures 8-10. It is particularly important that we investigate the response of the system to variation of A due to the fact that we could not parametrize A from the literature. It can be observed that as A . The feedback effects of buckling on growth, as represented by the term C 3 = C n 3 + C z 3 . Both C n 3 and C z 3 work to accentuate the magnitude of the buckle. The lines represent contours of constant R (in reference coordinates). The feedback owing to C n 3 works as follows. The cells positioned near the ellipses will receive more nutrient owing to the buckle, and are therefore more likely to proliferate. The increased proliferation of these cells leads to more compressive pressure that accentuates the buckle. Cells positioned near the squares receive less nutrient as a result of the buckle, and are therefore less likely to proliferate. Similarly, if we are incorporating anisotropic growth in our model, the feedback due to C z 3 works as follows. At the ellipses, the buckle relaxes the compressive stress tangential to the surface, which encourages more proliferation in this direction and accentuates the buckle.
increases from 0 to 20, the results often vary considerably. More precisely, many of the buckles change from being stable to unstable, and the magnitudes of the stable buckles often vary considerably. It seems likely that a major reason for this is that the radial displacement u 0 is positive at the time of buckling for low A, but as A increases the displacement u 0 at the time of buckling switches to being negative. This can be seen from the plot of the responsiveness of u 0 to A in figure 4.
We now discuss the effect of growth feedback, with reference to the plots in figure 8. Classical solid mechanics usually does not involve such feedback effects because the strain energy function is independent of the deformation. In the plots with growth feedback, we set k = 0 in equations (3.10) and (3.11), and in the plots without growth feedback, we set k = 1 in these equations. The growth feedback C 3 may be decomposed as C 3 = C n 3 + C z 3 , where C n 3 and C z 3 are given in (3.10) and the electronic supplementary material, equations (S4.18)-(S4.20). Here, C n 3 represents the growth feedback due to the effect of the buckle on the nutrient distribution, and C z 3 represents the growth feedback due to the effect of the buckle on the growth anisotropy. The growth feedback mechanism is explained in figure 11. In particular, in our results, C z 3 and C n 3 are both always negative. It can be seen, after an inspection of (3.11), that this means that both feedback terms accentuate the magnitude of the buckle. Furthermore, in light of the discussion in the electronic supplementary material, §S6 it means that the feedback terms do not destabilize the buckle.
In particular, when A = 0, the only feedback effect is due to C n 3 . The feedback due to C n 3 is almost always very small (i.e. |C n 3 | |C 2 |). This can be inferred from the fact that when A = 0 (so that the growth is isotropic and C z 3 is zero), the results with feedback are very similar to the results without feedback. However, as A increases, the growth becomes increasingly anisotropic and the relative magnitude of C z 3 increases. In fact, for large A, |C z 3 | |C 2 |, so that the equation (3.11) governing the magnitude of the buckle is dominated by the growth feedback term. This explains why, for large A, the results with feedback are very different from the results without feedback. Thus, it may be concluded that the feedback effect due to a perturbation in the nutrient distribution is negligible, but that the feedback effect resulting from the anisotropic growth is potentially very substantial.

Conclusions
The main thrust of this paper has been to develop a model that explains the buckling of a capillary through reference to growth in its immediate neighbourhood. As has been demonstrated by our 'infinite-geometry' results, the buckling can be explained purely as a consequence of the compressive stress arising from the growth gradient owing to the exponentially declining oxygen profile. We have found that planar buckling typically occurs on a time period of two to five times the cellular proliferation timescale, which is in good broad agreement with the literature. We have seen that feedback effects resulting from anisotropic growth can substantially accentuate the buckle, but that feedback effects resulting from a perturbation in the oxygen concentration are negligible. We expect that, generally, buckles of low azimuthal wavenumber are preferred by the system because they are usually energetically more favourable, and are less likely to be dissipated by the cellular phase. This is despite the fact that such buckles often go unstable last. Finally, it seems unlikely that our growth model can explain capillary tortuosity, because the required axial modes (with axial wavenumber 1) usually go unstable much later than the planar modes.