Nonlinear simulation of an elastic tumor-host interface

We develop a computational method for simulating the nonlinear dynamics of an elastic tumor-host interface. This work is motivated by the recent linear stability analysis of a two-phase tumor model with an elastic membrane interface in 2D. Unlike the classic tumor model with surface tension, the elastic interface condition is numerically challenging due to the 4th order derivative from the Helfrich bending energy. Here we are interested in exploring the nonlinear interface dynamics in a sharp interface framework. We consider a curvature dependent bending rigidity (curvature weakening) to investigate metastasis patterns such as chains or fingers that invade the host environment. We solve the nutrient field and the Stokes flow field using a spectrally accurate boundary integral method, and update the interface using a nonstiff semi-implicit approach. Numerical results suggest curvature weakening promotes the development of branching patterns instead of encapsulated morphologies in a long period of time. For non-weakened bending rigidity, we are able to find self-similar shrinking morphologies based on marginally stable value of the apoptosis rate.

area to volume ratio. In particular, regions where instability first occurs continue to grow at a faster rate than the rest of the tumor tissue. Such diffusional instability is induced by non-uniform cell proliferation and migration according to a heterogeneous distribution of nutrients. The tumor morphology is thus determined by the dominant nutrient levels where proliferation would be favored [Cristini et al.2005].
In the past several decades, mathematical models based on fluid mechanics were developed to understand the bio-mechanical properties of the tumor and its metastasis patterns. For example, the original model using Darcy's law (flow through a porous media) [Greenspan1976,Roose et al.2007] is composed of two parts. One is the concentration of a generic nutrient (e.g. oxygen or glucose) function σ satisfying a reaction diffusion equation σ t = D∆ σ − λ u σ , where D is the diffusion constant and λ u is the uptake rate; the other part is the internal pressure field p for tumor cell proliferation, which is related to cell velocity v by the Darcy's law v = −µ∇p, where µ is the cell mobility. The two is linked by the mass conservation of incompressible tumor cells ∇ · v = λ p (σ ), where λ p (σ ) = bσ − λ A is the cell proliferation rate. λ M = bσ ∞ and λ A are the rates of mitosis (cell birth) and apoptosis (cell death), respectively. This model is closed by introducing the Laplace-Young condition for internal pressure (p) ∂ Ω = γκ (derived from the surface energy of the tumor-host interface), far-field boundary condition for nutrient (σ ) ∂ Ω = σ ∞ , and the normal velocity of the moving interface V = −µ(∇p) ∂ Ω · n , where γ is the surface tension coefficient and κ is the mean curvature. The Darcy flow model, though simple in formulation, captures fundamental features of tumor mechanics and serves as a foundation for developing more sophisticated models, e.g. bifurcation behavior of the tumor growth by Stokes equation Hu2007a,Friedman andHu2007b].
In [Pham et al.2018], we introduced a two-phase Stokes model and treated tumor and its host as viscous fluids with different viscosity. The viscosities reflect the combined properties of cell and extracellular matrix mixtures. The tumor cell population is assumed to be homogeneous and cell proliferation produces a pressure field for tumor growth. Under the quasi-steady state assumption, this two-phase tumor model consists of a modified Helmholtz equation for nutrient diffusion in tumor and a Stokes equation for tumor dynamics. In [Turian et al.2018], we extended the two-phase Stokes model by introducing an elastic tumor-host interface governed by the Helfrich bending energy [Helfrich1973]. We derived a modified Laplace-Young condition of the stress jump across the interface for the Stokes equation using an energy variation approach, and performed a linear stability analysis to show how physical parameters such as viscosity, bending rigidity and apoptosis contribute to the morphological instability. Linear results suggest that increased bending rigidity versus mitosis rate contributes to a more stable morphological behavior, and there may exist fingering patterns for increasing tumor viscosity or apoptosis rate. Comparison with experimental data on glioblastoma spheroids shows good agreement, especially for tumors with high adhesion and low proliferation.
In this paper, we investigate the nonlinear dynamics of an elastic tumor-host interface. We introduce the bending rigidity coefficient as a function of local mean curvature, referred as curvature weakening model to describe broken intermolecular bonds and reduced stiffness of the interface [He et al.2012]. This is also motivated by results from recent studies that changes in stiffness of extracellular matrix may lead to increased mitosis and migration [Mason B.N.2012]. We reformulate the nutrient and Stokes equations as boundary integrals, and develop a sharp interface approach to explore the nonlinear instability of the interface. Note that in the sharp interface framework, we have to compute the 4th order derivative in the interface condition explicitly. After reformulation, the original two-dimension problem is reduced to one-dimensional curve integrals, which can be evaluated using spectrally accu-rate quadratures. To add more efficiency to the whole algorithm, a non-stiff interface updating scheme based on the small scale decomposition is implemented [Hou et al.1994]. Our numerical method is spectrally accurate in space and 2nd order accurate in time. Nonlinear simulations show that curvature weakening promotes the development of branching patterns and inhibits encapsulated morphologies. For non-weakened bending rigidity, there exist selfsimilar shrinking morphologies once the time dependent apoptosis rate (marginally stable value of the apoptosis rate) is applied. Though preliminary, the self-similar idea helps shed light on the strategy for morphological control, as a time dependent apoptosis might be enforced by a well-designed chemo-or radiotherapy.
This paper is organized as follows. In Sect. 2, we formulate the sharp interface model and non-dimensionalize the resulting PDE systems. In Sect. 3, we develop BIM formulation and present our numerical method including layer potential evaluations for boundary integrals and small-scale decomposition to remove stiffness. In Sect. 4, we show numerical results, and then we conclude with Sect. 5.

Mathematical model
We consider an avascular two-dimensional tumor as illustrated in Fig. 1. Let Ω 1 (t) be the tumor and Ω 2 (t) be the host tissue. The tumor-host interface Γ (t) is considered to be sharp and modeled as an elastic membrane. At the interface, a homogeneous elastic bending energy has been widely used to describe the interface dynamics either in a sharp or diffuse interface framework, see e.g. [Mikucki and Zhou2017, Adkins and Zhou2017, Gavish et al.2012, Du et al.2004,Du et al.2005,Wei2010,Dai and Promislow2013, Sohn et al.2012 among many others-i.e., the elastic bending energy, E H = 1 2 Γ (t) ν 0 κ 2 ds, where ν 0 is the constant bending rigidity coefficient, κ is the local mean curvature, and arc-length s parameterizes the interface. In general, one may consider a space dependent energy E W = 1 2 Γ (t) ν(κ)κ 2 ds, where the bending rigidity ν is given by a curvature weakening model [He et al.2012,Zhao et al.2016, where 0 ≤ C < 1 is the rigidity fraction and λ c is the characteristic radius beyond which ν(κ) decays significantly. When κ gets large, parameter ν(κ) approaches to its lower bound (1 −C)ν 0 . The largest rigidity limit ν 0 can be reached by setting C = 0. Performing domain variation in normal direction, i.e. compute variation δ E δΓ n := d dε E(x + εφ n), we obtain where the subscript s denotes a derivative with respect to the arclength parameter s and the prime notation is for a derivative with respect to κ. Nutrient field. Similar to the Darcy's model [Cristini et al.2003], here the nutrient field in Ω 1 (t) is governed by a reaction diffusion equation:

The two-phase Stokes tumor model
where D and λ u are the diffusion constant and uptake rate, respectively. For simplicity, we assume the nutrient concentration σ is constant in Ω 2 (t) and continuous across the interface Flow field. The mass conservation in Ω 1 (t) and Ω 2 (t) reads: where λ M and λ A are the mitosis and apoptosis rate, respectively. The Stokes equations in both domains are where T i = µ i ∇v i + (∇v i ) T + µ i (∇ · v i ) I − p i I are stress tensors for the interior tumor (i = 1) and exterior host (i = 2), p i are pressures, and parameters µ i and µ i are respectively the shear and bulk viscosity coefficients [Turian et al.2018, Pham et al.2018]. The stress jump condition across the interface is given by the bending energy variation in Eq. (2), [Tn] = δ E H δΓ n n or δ E W δΓ n n.
Note that the stress tensors take into account the rate of strain, dilatation and pressure. The normal velocity is simply V = v · n at the interface. Here we assume the cell velocity v is continuous across Γ (t), i.e.

Non-dimensionalization
Following [Turian et al.2018,Pham et al.2018, the dimensional variables are scaled by their characteristic values to yield the following non-dimensional parameters: , and σ ∞ are the characteristic diffusion length, time, and nutrient concentration scales, respectively. Also, P 1 = T 1 = µ 1 λ M . Since the tumor volume doubling time scale is typically much larger than the diffusion time scale (e.g. days vs. minutes), we assume λ M λ u , which leads to a quasi-steady reaction diffusion equation for nutrient field in the tumor tissue. Dropping all tildes, the nondimensional Stokes-flow system is given by: -In the tumor region Ω 1 (t), we have:
Conservation of tumor mass.
where A = λ A λ M represents the relative rate of cell apoptosis to mitosis. -In the host tissue region Ω 2 (t), we have a constant nutrient field: and divergence free condition for velocity and stress tensor: where T 2 = λ ∇v 2 + (∇v 2 ) T − p 2 I and λ = µ 2 µ 1 is the ratio between the exterior host and interior tumor viscosities.
-On the tumor-host interface Γ (t), we have no-jump boundary conditions for nutrient and velocity field: and a jump boundary condition for the stress tensor: where parameter S −1 = ν 0 µ 1 λ M L 3 represents the relative strength of bending rigidity.
Ideally, one would like to rewrite the Stokes system to the standard one, in which the velocity fields are divergence free in both the tumor and host regions. This is helpful in the design of numerical methods. To do this, we redefine the tumor cell velocity in Ω 1 (t) as where d = 2 is the spatial dimension. Using the modified Helmholtz equation (9) and the identity ∇ · x = d, equation (11) becomes divergence free: Thus the PDE system in Ω 1 (t) becomes: wherep 1 is the renamed interior Consequently, the boundary conditions can be rewritten as: where T u 1 = ∇u 1 +(∇u 1 ) T −p 1 I, f (κ) = 1 2 κ 3 +κ ss for E H and f (κ) = 1 2 ν + 2ν κ + ν κ ss + 1 2 ν κ 2 + 3ν κ + 3ν κ 2 s + 1 2 ν κ + 1 2 ν κ 3 for E W . The reformulation requires an explicit evaluation of ∇∇σ n, which can be expressed in terms of the normal derivative of σ by Eq. (9)as where s is the arclength representation of the tumor-host interface. Since the exterior velocity field is already divergence-free, it is unnecessary to reformulate the exterior problem. Note that the reformulated velocity field becomes discontinuous across the interface. Notice that although in Eqs. (19), (20) both the velocity and the nutrient field are governed by time-independent PDEs, the boundary itself and the boundary conditions in Eq. (21) are time-dependent, which results in a moving boundary problem.

Review of linear stability analysis.
Though a linear stability analysis could be performed using a weakened bending rigidity, the calculation is very tedious and cumbersome. For brevity, we focus on a constant bending energy case [Turian et al.2018] and mainly use our numerical solvers to study the curvature weakening model. For a circular tumor spheroid of radius R(t), the interface evolves as: where I 0 (R) and I 1 (R) are the modified Bessel functions of the first kind with indices 0 and 1, respectively. Figure 2 shows the rate of change of a tumor spheroid with respect to the radius R for different A . R s (A ) is the linear steady radius associated with A satisfying dR/dt = 0.
Large A causes more cell death and therefore limits the size of the tumor spheroid, whereas A = 0 indicates an unbounded growth. For a slightly perturbed circular interface, r(α,t) = R(t) + εδ (t) cos(lα), where α is the polar angle, δ (t) is the time-dependent perturbation, ε 1, and integer l ≥ 2 is the perturbation mode. The shape perturbation for Stokes-flow with a uniform bending energy evolves as where the shape factor δ R measures the deviation of the tumor shape from a circle of varying radius, thus describing tumor morphological stability. Note that the normal contribution of the term ∇∇σ in the stress jump is focused here. The shape perturbation depends on A , S −1 , l and λ . Observe that the right hand side of Eq. (24) increases with increasing A (high cell death) and decreases with increasing S −1 (high membrane rigidity), implying that A promotes shape instability while S −1 stabilizes it. The morphological stability is mainly determined by the competition between these two important parameters. The parameter λ may promote or reduce instability, depending on the values of l and the radius of the tumor R. In the result section, we will do a parameter study on the critical value of the stiffness S −1 and examine the full nonlinear dynamics.

Numerical method
We use the boundary integral methods to solve (1) the nutrient field in tumor domain; (2) the 2D Stokes equation for the fluid velocity field in both domains. We then update the position of the interface Γ (t) by a nonstiff 2nd order multistep method. The algorithm presented below is an extension of the approach developed in [Pham et al.2018] for interfacial flows with surface tension. For completeness, we outline the main ideas here. A rigorous convergence and error analysis of the boundary integral method for the tumor problem can be found in [Hao et al.2018].

Boundary integral formulation
The nutrient field. Consider a Green's function for the modified Helmholtz equation in where G σ = G σ (x, x ), x ∈ Ω 1 (t) is the source point, x is the field point, and δ x (x, x ) is the Dirac delta function. Thus the fundamental solution for Eq. (9) is the modified Bessel function of the second kind where r = |x − x |. Using potential theory [Kress2013], we define a single-layer potential for the modified Helmholtz equation: where ζ is the layer potential. The nutrient σ can thus be written as a double-layer potential: where n is the unit outward normal to Γ (t). By the uniform nutrient condition (21), we may repose Eq. (9) as a second-kind Fredholm integral equation with an unknown density ζ on Γ (t): The term ∂ σ ∂ n (normal derivative of double-layer potential) can be computed using The Stokes-flow field. Following [Pham et al.2018], let G be the Stokeslet and T be the tensor stresslet, then the single layer potential at the interface Γ is and the double layer potential at Γ is where P.V. indicates the principle value integral. Assuming that the flow vanishes at the farfield, the boundary integral representation of the velocity v 2 approaching the interface from the exterior domain and the velocity u 1 approaching the interface from the interior domain Ω 1 is where f i denote the interior and exterior normal stress at the interface λ −1 T 2 n = f 2 and T u 1 n = f u 1 ). Since we do not know T i individually, we rewrite Eqs. (33) and (34) in terms of T 2 − T u 1 and make use of Eq. (21). To start, we multiply Eq. (33) by λ and add it to Eq.(34) to get where the term T 2 − T u 1 is explicitly given in Eq. (21). At the interface, from Eqs. (14), (17), the interior and exterior velocities are related by Putting Eqs. (35) and (36) where the force term Note that in 2D, v 2 = (v 1 , v 2 ) , n = (n 1 , n 2 ) , and F = (F 1 , F 2 ) . Using the formulas of the double and single layer potentials [Pozrikidis1992], equation (37) can be explicitly rewritten as Hence, the explicit forms of the single layer and double layer potentials are and s (α ))) and n i = n i (x(s(α))). In Eq. (41), the only singularity in the integrand comes from the logarithmic kernel. This can be analyzed in the following subsection.

The evaluation of the boundary integrals [Li and Li2011]
With the integral formulation above, we assume the interface Γ is analytic and given by x(α,t) = (x(α,t), y(α,t) : 0 ≤ α ≤ 2π , where x is 2π-periodic in the parametrization α. The unit tangent and normal (outward) vectors can be calculated as s = (x α , y α )/s α , n = (y α , −x α )/s α , where the local variation of the arclength s α = x 2 α + y 2 α . Subscripts refer to partial differentiation. We track the interface Γ by introducing N marker points to discretize the planar curves, parametrized by α j = jh, h = 2π N , where N is a power of 2.
Computation of the single-layer potential.
In Eqs. (27) and (41), the single-layer potential type integrals contain the Green functions with the logarithmic singularity at r = 0. They can be rewritten as the following form under the parametrization α: where Φ are the Green functions G or G σ , φ is the layer density η or ζ . We decompose the Green functions as: where I 0 is a modified Bessel function of the first kind, r = |x(α) − x (α )|. The square bracket on the right-hand side of Eqs. (44) and (45) has removable singularity at α = α , denotes a smooth function that vanishes as α → α , and since K 0 has the expansion As a result, for an analytic and 2π-periodic function f (α, α ), a standard trapezoidal rule or alternating point rule can be implemented to evaluate the integral and achieve spectral accuracy. The first term on the right-hand side of Eqs. (44) and (45) is still singular and can be evaluated using the following spectrally accurate quadrature [Kress1995]: where m = N 2 , α i = πi m for i = 0, 1, ..., 2m − 1, and weight coefficients The derivative d dα in Eq. (43) is approximated using Fast-Fourier-Transform spectral derivatives thus maintaining spectral accuracy.
Computation of the double-layer potential. In Eq. (28), the double-layer potential type integrals contain the Green functions with logarithmic singularity at r = 0. It can be rewritten as the following form under the parametrization α: where Φ stands for the Green function G σ and φ is the layer density ζ . Since ∂ G σ ∂ n has logarithmic singularity, we decompose it as below: where g 1 (α, α ) and g 2 (α, α ) are analytic and 2π-periodic functions with We have used the fact d dr K 0 (r) = −K 1 (r). Since K 1 has the expansion , the square bracket on the right-hand side of Eq. (52) also has a removable singularity at α = α , thus the integral involving g 2 (α, α ) can be evaluated by a standard trapezoidal rule or alternating point rule. Note that The first term on the right-hand side of Eq. (52) is still singular and can be evaluated through the quadrature given in Eqs. (47) and (48).
To summarize, using Nyström discretization with the Kress quadrature rule [Hao et al.2014] described above, we discretize the boundary integral equations for the nutrient and Stokes fields into two dense linear systems with unknowns as the layer densities ζ and η on Γ (t), which can be solved using an iterative solver, e.g., GMRES [Saad and Schultz1986].

The evolution of the interface
As indicated in [Hou et al.1994], the curvature driven motion introduces high-order derivatives, both non-local and non-linear, into the dynamics through the Laplace-Young condition at the interface. Explicit time integration methods suffer from severe numerical stability constraints and implicit methods are difficult to apply since the stiffness enters non-linearly. Hou et al. resolves these difficulties by adopting the θ − L formulation and the small-scale decomposition (SSD) which we will follow in this paper.
θ − L formulation. This description makes the application of an implicit method straightforward and may circumvent the problem of point clustering. Consider a point x(α,t) = (x(α,t), y(α,t)) ∈ Γ (t).
Denote the normal and tangent velocity by V (α,t) = u · n and T (α,t) = u · s respectively, where u = x t = V n + T s describes the motion of Γ (t). The tangent angle that the planar curve Γ (t) forms with the horizontal x-axis, called θ , satisfies θ = tan −1 y α x α . The unit tangent and normal vectors become s = (cos θ , sin θ ) and n = (sin θ , − cos θ ). The length of one period of the curve is L(t) = 2π 0 s α dα. Differentiating these two equations of θ , s α in time, we obtain the following evolution equations: where the curvature is evaluated by κ = θ s = θ α s α . Instead of using the (x, y) coordinates, we are able to repose the equation of motion in terms of dynamical variables (L, θ ). To gain more efficiency and accuracy, one may choose a tangent velocity T (independent of the morphology of the interface) such that the marker points are equally spaced in arclength to prevent point clustering: It follows that s α is independent of α and thus is everywhere equal to its mean: The procedure for obtaining the initial equal arclength parametrization is presented in Appendix B of [Baker and Shelley1990]. The idea is to solve the nonlinear equation for α j using Newton's method and evaluate the equal arclength marker points x(α j ) by interpolation in Fourier space. We may recover the interface by simply integrating Small scale decomposition (SSD).
The idea of small scale decomposition (SSD) is to extract the dominant part of the equations at small spatial scales [Hou et al.1994]. To remove the stiffness, we use SSD in our problem and develop an explicit, non-stiff time integration algorithm. Through the analysis of the single-layer and double-layer terms, the only singularity in the integrands comes from the logarithmic kernel. Following [Hou et al.1994] and noticing the curvature terms in the stress-jump condition in Eq. (21) and Eq. (39), one can show that at small spatial scales [Sohn et al.2012], where H (ξ ) = 1 2π 2π 0 ξ cot α−α 2 dα is the Hilbert transform for a 2π-periodic function ξ . We rewrite Eq. (54), where the Hilbert transform term is the dominating high-order term at small spatial scales, , where ∆t and h are the time-step and spatial grid size, respectively. The efficiency has been demonstrated numerically in the seminal work [Hou et al.1994] and later in [Li et al.2007,Zhao et al.2017] for a Hele-Shaw problem. For the tumor growth problem, the semi-implicit time-stepping scheme (see Eq. (60)) requires ∆t = O(h) instead of explicit schemes which would require ∆t = O(h 3 ). In section 4, we show numerical examples using N = 2048. In the simulation, we could use ∆t as ∆t = 1.0 × 10 −2 for stability instead of ∆t < 10 −6 for an explicit scheme.

Semi-implicit time-stepping scheme
Taking the Fourier transform of Eq. (61), we get In Fourier space, we solve Eq. (62) using a second order accurate linear propagator method in the Adams-Bashforth form and then apply the inverse Fourier transform to recover θ . Specifically, we discretize Eq. (62) aŝ where the superscript n denotes the numerical solutions at t = t n and the integrating factor . (64) Note that by setting the integrating factors in Eq. (63) to 1, we recover the classical Adams-Bashforth explicit time-stepping method. The integrating factors in Eq. (63) can be evaluated simply using trapezoidal rule, To compute the arclength s α , equation (55) is discretized using the explicit 2nd-order Adams-Bashforth method, where M is calculated using M = 1 2π 2π 0 V (α,t)θ α dα. Note that the second order linear propagator and Adams-Bashforth methods are multistep methods and require two previous time steps. The first time step is realized using an explicit Euler method for s 1 α and a first order linear propagator of a similar form forθ 1 . To reconstruct the tumor-host interface (x(α,t n+1 ), y(α,t n+1 )) from the updated θ n+1 (α) and s n+1 α , we first update a reference point (x(0,t n+1 ), y(0,t n+1 )) using a second-order explicit Adams-Bashforth method to discretize the equation of motion x t = Vn (with the tangential part dropped since it does not change the morphology) (x(0,t n+1 ), y(0,t n+1 )) = (x(0,t n ), y(0,t n )) + ∆t 2 (3V (0,t n )n(0,t n ) −V (0,t n−1 )n(0,t n−1 )) .

Results
We have computed a number of different cases which illustrate and expand upon the linear stability analysis in [Turian et al.2018]. First, we examine stability of a perturbed interface through the full nonlinear simulations. We then take into account the curvature weakening of the bending rigidity and examine its effects on the pattern formation. Finally, we examine the existence of possible self-similar solution.
The correctness of implementations of boundary integral methods for both Stokes flow and the modified Helmholtz equation was checked in a number of ways. This includes growing/shrinking of a circular interface to the steady radius predicted by Eq. (23). Agreement of the nonlinear evolution of perturbations was also verified with linear solution for a slightly perturbed circular interface in Eq. (24). This is assessed by comparing the corresponding linear and nonlinear shape factors. The linear shape factor is calculated by solving Eqs. (23) and (24), and the nonlinear shape factor is calculated numerically using where x j denotes the discrete points that describe the tumor/host interface and R denotes the effective radius of the tumor, which is the radius of a circle with the same area as the tumor.
Recall that S −1 is proportional to membrane rigidity. In Fig. 3 [a], we illustrate this behavior by plotting S −1 with A = 0.5, mode l = 3 for various viscosity ratios λ = 0.5, 1.5, 2.5 against R. When A = 0.5, the steady radius is R s ≈ 3.326. We consider the dynamics of tumors that may grow or shrink depending upon their initial radius. We take two membrane rigidity parameters S −1 = 0.001 and S −1 = 2, two initial tumor radii R(0) = 1.988 and R(0) = 4.5, and vary the viscosity ratio λ .
Specifically, we consider in Fig. 3 [a] evolution from the points P 1 (1.988, 0.001), P 2 (4.5, 0.001), Q 1 (1.988, 2), Q 2 (4.5, 2) where the first coordinate represents the initial tumor radius R(0) and the second represents the membrane rigidity parameter S −1 . When A = 0.5, linear theory predicts that a circular interface will evolve to its stationary radius R s ≈ 3.326, while the stability of a perturbed interface depends on S −1 . As seen from Fig. 3 [a], linear theory predicts that starting from the point P 1 (1.988, 0.001), P 2 (4.5, 0.001), where membrane rigidity is low, the 3-mode perturbation will be unstable as the tumor grows or shrinks to the stationary radius R s . We point out that the effective tumor radius will actually grow/shrink to the stationary radius R s (A = 0.5) ≈ 3.326 for a while and turn out to be larger than this predicted size due to morphological instability. However, starting from the point Q 1 (1.988, 2), Q 2 (4.5, 2) where the membrane rigidity is higher (S −1 = 2), the simulation shows tumor will grow or shrink to R ≈ 3.326 and is stable.
In Fig. 3 [b] and [c], we show the nonlinear evolution of tumors starting from points P 1 and Q 1 , respectively. In each case, the initial tumor/host interface is a 3-mode perturbation of a circle given by r = 1.988 + 0.05 cos(3α). Similar computations for P 2 and Q 2 are shown in Fig. 3 [d] and [e] with initial shape r = 4.5 + 0.05 cos(3α). The viscosity ratio is varied from λ = 0.5, 1.5, 2.5 as labeled in the plots. The nonlinear shape factors are plotted as functions of time, and the corresponding tumor morphology at the final time t f of each simulation are shown as insets. We also plot the linear solution for λ = 0.5 (dashed line) in Fig. 3 [b] and [d]. Note that the calculations for Fig. 3 [b] and [d] stop because higher numerical resolution is needed to resolve high curvature regions.
When the membrane rigidity is small (S −1 = 0.001), Fig. 3 [b] show that the shape perturbation starting from the point P 1 grows rapidly, especially for small viscosity ratio, indicating an unstable growth. This is consistent with the predictions of linear stability theory ( Fig. 3[a]). On the other hand, when the membrane rigidity is increased to S −1 = 2, nonlinear simulations from point Q 1 converge to a final circular morphology for all three cases and the tumor grows stably to its diffusion limited size R s ≈ 3.326, as shown in Fig.  3[c]. Simulations of points P 2 and Q 2 in Fig. 3 [d] and [e] show similar unstable and stable shrinking behavior, respectively.
In Fig. 3 [f], using the point P 1 , we compare the shape perturbation between the nonweakening bending model (constant bending stiffness plotted using solid lines) and curvature weakening bending model (plotted using dashed lines). Here we set λ c = 1.25,C = 0.95 in Eq. (1). For small viscosity ratio λ = 1.5, the curvature weakening effect dramatically slows down the process of unstable growth and leads the interface morphology to branching patterns. However, as λ increases, the weakening effect is reduced and we get the usual encapsulated morphology. These results highlight the level of complexity and sensitivity of interface dynamics in fluid due to inhomogeneous elasticity.

Growth in nonlinear regime (A = 0.7)
We next increase the apoptosis rate A from 0.5 to 0.7 and focus only on the viscosity ratio λ = 1.5 case. Beyond the parameter regime of linear prediction, in Fig. 4 [a] the interface evolves far away from the steady radius R s (A = 0.7) = 1.988. Notice that in this case based on its initial 3-mode perturbation, i.e. r = 1.988+0.05 cos(3α), the interface develops further splitting and fingering patterns due to the curvature weakening effect as shown in Fig.  4 [c]. The evolving morphology changes from a compact shape with inward splitting shown in Fig. 4 [b] to a fingering pattern with a tendency of outward splitting as shown in Fig. 4 [c]. Here we use curvature weakening parameters λ c = 1.25, C = 0.95.
We also compute the evolution of a complex initial shape in Fig. 5, where r = 1 + 0.05 1.988 cos(2α)+ 0.1 1.988 cos(3α)+ 0.08 1.988 sin(4α)+ 0.12 1.988 cos(5α). In Fig. 5 [a], tumor evolution with complex initial shape grows unstably with appearance of long and slim zigzags inside the tumor; while for the one with curvature weakening effect (λ c = 1.25, C = 0.95) in Fig. 5 [b] such pattern disappears, and it takes much longer time to reach the size in Fig.  5 [a].
of a self-similar sequence are plotted in Fig. 6 [b] and [c]. Although preliminary, the selfsimilar idea helps shed light on the strategy for morphological control, as a time dependent apoptosis might be enforced by a well-designed chemo-or radiotherapy.

Conclusion
In this paper, we performed nonlinear simulations of a 2D, non-circular tumor with isotropic or curvature weakened bending rigidity growing in a host tissue. The interior tumor and exterior host were modeled by the Stokes flow, and the tumor-host interface was modeled by an elastic membrane governed by the Helfrich bending energy. Using boundary integral formulations of the Stokes flow and the nutrient field, we developed a spectrally accurate sharp interface method. We then investigated the nonlinear dynamics of the tumor-host interface.
The linear stability analysis suggests that an increase in bending rigidity contributes to an increase in morphological stability for an isotropic bending rigidity. Nonlinear simulation confirms this and moreover, curvature weakening bending helps improve the stability by slowing down the growth of shape perturbations and promotes branching or tip-splitting fingering patterns rather than encapsulated morphologies for small viscosity ratio, In fact, not only for the Stokes model, our recent preliminary results using the Darcy's model suggest more pronounced fingering patterns if a curvature weakening bending is implemented, as shown in Fig. 7. We can see the self-branching morphology is enhanced with curvature bending energy as reported in [Zhao et al.2016] for a Hele-Shaw interface. It is also observed that an increase in the apoptosis leads to an overall increase in shape instabilities.
In experiments, thermal and mechanical stresses have been found to be important in regulating cell fates and motility, proliferation and apoptosis rates. In future work, we will consider these effects. We will also consider adding a stochastic component to the current model. Although our 2D results are expected to hold qualitatively in three dimensions as suggested by the linear stability analysis (at least for Darcy's model), we would like to perform full 3D simulations to confirm this.
We have chosen the initial shape as a perturbed circle since that that in vitro tumor grows nearly spherical at early times. It is reasonable to assume that, in vivo, tumor at its initial stage of avascular growth is nearly spherical. As the tumor continues to grow, morphological instabilities come in and the tumor starts to develop protruding shapes. In this work we focus on the morphological instability of the interface and our numerical scheme can handle even complex tumor morphologies as indicated in Fig. 4[a] and Fig. 5[a] with spectral accuracy in space. For studies involving simulation beyond topological changes, we plan to use phasefield or level-set formulation, which is our future work. . Apoptosis A is time-varying such that the shape factor δ R = 0 in linear regime. Initial shape: r = 2 + 0.2 cos 3α for [b] and r = 3.5 + 0.35 cos 3α for [c]. [a] [b] Fig. 7 Growth with Darcy-flow model with [a] isotropic bending energy and [b] curvature weakening bending (stiffness fraction C = 1, characteristic length λ c = 1); the parameters are set as: mesh points N = 2048, time step dt = 0.01, bending rigidity S −1 = 0.15, apoptosis A = 0.7, initial shape: r = 2 + 0.05 cos(2α).