Fingering instabilities in tissue invasion: an active fluid model

Metastatic tumours often invade healthy neighbouring tissues by forming multicellular finger-like protrusions emerging from the cancer mass. To understand the mechanical context behind this phenomenon, we here develop a minimalist fluid model of a self-propelled, growing biological tissue. The theory involves only four mechanical parameters and remains analytically trackable in various settings. As an application of the model, we study the evolution of a two-dimensional circular droplet made of our active and expanding fluid, and embedded in a passive non-growing tissue. This system could be used to model the evolution of a carcinoma in an epithelial layer. We find that our description can explain the propensity of tumour tissues to fingering instabilities, as conditioned by the magnitude of active traction and the growth kinetics. We are also able to derive predictions for the tumour size at the onset of metastasis, and for the number of subsequent invasive fingers. Our active fluid model may help describe a wider range of biological processes, including wound healing and developmental patterning.

What causes the formation of such fingers? Studying cancer has traditionally focused on a large number of biological (especially genetic and biochemical) cues [13]. Yet essentially operate by collectively affecting a smaller number of physical properties of the tissues and environments involved [1]. How these physical alterations can, in turn, lead to fingering has been investigated by several models. Various causal mechanisms, differing in their assumptions on the mechanical properties adopted for the tissues, have been proposed [14 -20]. For example, reactiondiffusion models of nutrient-limited growth have been used [1,21] because the accessibility of diffusing chemicals is necessary for tumour growth. But, other models have treated the fingers' emergence as a consequence of a mechanical instability [14,15,17,22]. Support for the latter approach was provided by experiments conducted to prevent biochemical signalling, but in which fingering occurred nevertheless [23]. Mechanical processes proposed by these models include: fracturing of an elastic surrounding medium (extracellular matrix, healthy tissue) upon pressure from a growing solid inclusion (tumour), leading to subsequent infiltration of the cracks by malignant cells [22]; branching instabilities caused by the tumour's confinement pressure, analogous to the corona splash of a drop impacting a solid surface [24]; mechanical frustrations between an outer, proliferating ring of a growing tumour and its necrotic core [15]; buckling due to swelling of a spatially restricted gel-like tissue [17]; instability resulting from the interplay between spatially non-uniform cell division/death rates and shear in viscous tissues [14]; or a pulling mechanism by a subpopulation of leader cells at the tumour's edges [10,25].
Many of the above mechanisms have been successfully described using continuum, analytically solvable models [14 -17], thus providing deep insights into the physical context involved in tissue fingering. However, these continuum models of tumours have, in most cases, omitted one fundamental component of live tissues: cells actively apply forces (via the conversion of chemical energy into mechanical energy) to their surroundings [10,26 -30]. Yet, the onset of fingering at tissue boundaries is often correlated with the dense presence of these so-called active forces and the occurrence of the resulting self-propelled motion [10,29]. In discrete models, the central role of these forces in triggering tissue fingering has been verified by simulations within several different frameworks [18,25]. But, to the best of our knowledge, only few continuum models of fingering have studied the consequences of self-propelling forces: the effect of an active rim at the boundary of the tissue has been studied by Mark et al. [31] and by Nagilla et al. [32], while the role of active forces in the tissue bulk has been discussed in wound healing models by Zimmermann et al. [33] and by Nesbitt et al. [34]. However, the analyses in the two latter works are limited to a non-dividing, rectangular tissue, either in a static state or somehow pushed on one side by a rigid barrier.
We here investigate the role of tissue bulk activity for the emergence of fingers in more general situations. We first construct a continuum mechanical model of an active and growing tissue, supported by experimental evidence, and that is analytically solvable and that involves only four physical parameters: friction, activity, growth and surface tension. We next investigate the role of activity in promoting fingering. Provided with experimentally derived estimates of the physical parameters, our model notably produces realistic predictions for the number and evolution of the fingers.

Assumptions and equations
In our model of an active and growing tissue, the evolutions of the pressure p(r, t) and velocity v(r, t) fields at a position r and time t are governed by the following force balance and mass conservation equations, and respectively, with r the nabla operator, and where a and b are positive parameters, specifying the strength of the interaction between the tissue and a substrate: a describes the magnitude of the active traction, while b describes the magnitude of the effective passive friction ( proportional to the tissue viscosity). In equation (2.1b), k is the net rate of growth (we are here interested in regimes in which it is also positive) of an incompressible tissue, undergoing cell division (or individual cell growth). The a-term in equation (2.1a), which accounts for the tissue activity, is discussed in detail in §2.2. The evaluation of the various parameters is examined in §2.3. Ignoring the a-term, equation (2.1a) reduces to Darcy's Law rp ¼ Àbv (originally used to describe viscous flows in porous materials and Hele -Shaw apparatus [35]), which has been widely employed to model the passive behaviour of tissues [1,36 -38]. Darcy's Law notably assumes a viscous and quasi-two-dimensional dynamics for the deformations of a tissue layer, by considering that the effects of friction against a substrate are much stronger than those of viscous shear within the plane of the layer.
Using two-dimensional models is experimentally justified by the large prevalence of in vitro tissue culture monolayers, and also because many in vivo soft tissues, including epithelium in which carcinomas develop, tend to spontaneously form quasi-two-dimensional monolayers [10,29,39]. Consequently, two-dimensional descriptions are often employed in tissue mechanics models [14,17,25].
The mechanical properties of live tissues at short timescales, up to the order of minutes, are generally dominated by an elastic constitutive behaviour. At longer timescales, however, a viscous description is better suited [30]. The crossover between the two regimes is probably related to the turnover rates of intercellular adherent junctions [40,41]. Hence, epithelial tissues become fluidized by a reduction in the number of adherent junctions, and a concomitant increase in the magnitude of active traction when becoming malignant. This well-known 'melting' process is often referred to as the 'epithelial to mesenchymal transition' [3,5,42]. Moreover, cell division and apoptosis are also known to further contribute to the fluidization of tissues [43]. As we here model the behaviour of the tissue at timescales on which it experiences substantial growth (that is, on the order of several hours at least [4]), the viscous constitutive behaviour implied by Darcy's Law is justified. Many existing continuum models of epithelial tissues indeed make the same assumption [14,30,38,44]. Note that, by writing equations (2.1), we further assume that inertial terms are negligible on these timescales.

Tissue activity
The second term on the right-hand side of equation (2.1a) accounts for cells actively propelling themselves by exerting traction against the substrate. It will subsequently be referred to as the active term, and a specifies its strength.
We consider here that the direction of the net local active force acting on the tissue layer from the substrate is aligned with the direction of the local flow velocity. This assumption was made in previous studies modelling active tissues [18,45]. It is a consequence of cells attempting to maintain their direction of motility, as illustrated in figure 1a, and also manifested by the persistent Brownian motion of individual cells in vitro [47]. On a subcellular level, it probably results from the friction destabilizing lamellipodia [11,48] that are not aligned with the cell's velocity [18]. It has further been shown to be consistent with experimental data [18,49,50], although this directionality is not universal [51] and other models have been proposed where the direction of the active force is treated as an independent internal variable, coupled to stress and velocity fields [44,52]. Equation (2.1a) also assumes that the active traction does not depend on the magnitude of the velocity. This assumption has been made in several numerical models of motile cells [18,53,54], and enables a distinct analysis of the role played by activity.
In the classic theory derived by Toner & Tu [55,56], often used to model active fluids [33,34,57], the net force per unit volume acting on the active fluid from the substrate, in a spatially uniform flow, follows F T (v) ¼ a T v 2 b T jvj 2 v (ignoring here additional inertial and gradient terms; with a T and b T positive parameters). In comparison, our model, equation (2.1a), gives this force the expression F(v) ¼ av/jvj 2 bv. In both cases, the fluid has a 'preferred' spontaneous magnitude of velocity jv s j ¼ a/b (¼ (a T /b T ) 1/2 in the Toner-Tu model), which it would select when moving in unbounded space without being driven by an external pressure gradient or growth. For magnitudes of velocity lower than jv s j, the fluid would be driven to move faster by the a-term, while above it, it would be slowed down by the friction (the b-term). This fluid's constitutive behaviour may be described in terms of a Landau-type 'velocity potential' [46], shown in figure 1b. Although our model has a discontinuity in the direction of F(v) at v ¼ 0 which does not exist in F T (v), we have verified that this singularity does not significantly affect our subsequent results.
The Toner -Tu model, however, assumes that the friction force grows as jvj 3 and that the active force varies linearly with jvj. Both assumptions are unrealistic when describing biological tissues. Our model, rsos.royalsocietypublishing.org R. Soc. open sci. 5: 181579 however, retains the physical interpretation of b as the friction coefficient of Darcy's Law, and of a as the magnitude of the active force of the tissue against the substrate ( per unit volume). Our approach also enables a direct comparison with classical results for viscous fingering [58], readily obtained from our model by taking the limit a ! 0.

Estimation of parameters
Based on in vivo microscopy observations [4] of the time necessary for doubling a carcinoma's size, which is on the order of a few hours, we estimate that the growth rate k is about 10 24 s 21 .
The passive friction b can be estimated based on in vitro force measurements of epithelial tissues against substrates [29,30,59] (admittedly, inferences about in vivo systems from these in vitro experiments are arguable). Following Pompe et al. [59], we assume that friction with the substrate is primarily the consequence of cell-substrate ligands, numbering 200-300 per cell, each of which exerts a force of about 10 212 N. We thus estimate the total friction force per cellular volume to be about 10 6 N m 23 for a 10 25 m cell size. From in vivo microscopy of micrometastasis growth [4], the typical velocity v of the cells falls within 10 210 to 10 29 m s 21 . Hence, dividing the volumic friction force by this velocity provides an estimate b 10 15 -10 16 Pa s m 22 .
There is no lower limit on a, as epithelial cells may not exert any active force against the substrate. The upper limit can be estimated on the basis of force tracking microscopy applied to spreading epithelial monolayers in vitro [29,30]. It is observed that traction forces are actively exerted through the monolayers and peak at their edges, giving rise to a gradient of the stress tensor's diagonal terms, which is up to 10 7 Pa m 21 in the study by Trepat et al. [29], and 10 8 Pa m 21 in the work of Blanch-Mercader et al. [30]. Balancing a with this typical stress gradient, one can place an upper estimate on a at approximately 10 8 Pa m 21 . In agreement, the traction exerted by single fibroblasts has been reported as up to 10 27 -10 25 N per cell [60], which would correspond to a 10 8 -10 10 Pa m 21 when dividing by the cell's volume.
The effective surface tension of the tissue g will also play a role in our further considerations. Its magnitude depends on the strength of intercellular adhesion and behaviour of cortical actin networks [61,62]. The surface tension was evaluated indirectly by Foty et al. [63] by measuring the energetic penalty of compression of embryonic multicellular spheroids, revealing values on the order of 3-9 mPa m. We thus presume that g 10 23 -10 22 Pa m is a realistic range for our system. A summary of the estimates for the physical parameters is presented in table 1.

Model system
The tumour is modelled as an initially circular, two-dimensional droplet of a growing active fluid described by equations (2.1), with an unperturbed, time-dependent radius r 0 (t) (figure 2). The surrounding healthy tissue is modelled as a passive, non-dividing fluid, whose pressure p 0 (r, t) and velocity v 0 (r, t) fields follow the equations: where b 0 is the friction parameter (analogous to b in the active fluid). In writing equation (2.2b), we effectively assume that growth in the passive fluid can be neglected on the timescale of metastasis initiation. We assume that the activity a and the growth rate k are constant (independent of r and t) through the active tissue. Constant magnitude of the active force has been assumed in models of active matter before [64,65]. While not correct in all situations [30,66], it is a convenient assumption to evaluate the influence of its magnitude in fingering. Extensions to non-uniform and time-dependent behaviours of these parameters are readily possible, and we investigate a case of evolving growth rate in §3.4 (see also appendix B.2). Similarly, we assume that b and b 0 are uniform within their respective regions.
We study the system in polar coordinates r ¼ (r, u) and write vector fields' components in this system with appropriate subscripts, such as v ¼ (v r , v u ). The perturbed interface between the active and passive fluids, described by the line r(t, u), must satisfy two boundary conditions. First, the continuity of the radial components of velocities is expressed as follows:  Second, the pressure difference across the interface separating the two tissues must equal the Laplace pressure, pj r¼r À p 0 j r¼r ¼ Àg with g the surface tension, and the fraction being the expression of the local interfacial curvature in polar coordinates [67]. Note that our description does not account for an interfacial bending force [31,68]. We indeed assume that the corresponding energetic cost, which involves a higher order dependency in interfacial curvature when compared with surface tension [31], is negligible at the length scales considered for the linear stability analysis presented in §3.
Unless stated otherwise, we use the dimensionless variables defined as follows. Distances are re-scaled by the characteristic length ' ¼ ð 2g bk Þ 1=3 , which can be interpreted as a capillary length at which growth balances interfacial tension (on the order of 10 mm, based on the estimates of table 1). Times are rescaled by k 21 , and we further define f ¼ b 0 /b the relative viscosity of the displaced tissue compared to the active growing droplet. We introduce a reference activity a Ã ¼ b'k 10 7 Pa m À1 to make the active traction a dimensionless, a=a Ã ! a. Pressures and velocities are made dimensionless by p Ã = b' 2 k 10 2 Pa and v Ã = 'k 10 À9 m s À1 , respectively. We will use the same letters for the dimensionless versions of the variables as for their dimensional counterparts.
The droplet of the active tissue grows due to a positive k, as required by equation (2.1b), and the passive fluid is displaced by it. As long as the interface between the two tissues remains circular (with the unperturbed radius r 0 ), hydrodynamic fields in both regions remain symmetric under rotations and are given by as obtained by solving equations (2.1)-(2.4) and using the dimensionless quantities defined above.

Linear stability analysis
We investigate under which conditions the active, circular droplet of radius r 0 would start to form fingerlike protrusions at its edge, while undergoing uniform growth. To do this, we perform a linear stability analysis around the circular solution given in equations (2.5), by investigating infinitesimal interfacial perturbations of the form r ¼ r 0 þ dr, with for an integer n corresponding to the mode of the periodic perturbations, and where f n (t) is a function describing its time evolution (with lim t!0 f n (t) ¼ 1 for all n). An analogous ansatz of periodicity in u is made for the perturbations of hydrodynamic fields in both fluids, (dp, dp 0 , dv, dv 0 ) / f n (t) e inu , around the solution given by equations (2.5). Note that our stability analysis concerns small perturbations around solutions that are themselves time-dependent. Applying the evolution equations, equations (2.1)-(2.4), to these perturbed fields provides an expression for the n-mode's rate of growth defined by Positive values of s n (r 0 ) correspond to unstable, growing modes n, which may become the basis for the formation of fingers. We derive in appendix A the following expression for s n (r 0 ): where L n (x) is the function with n j the binomial coefficient 'n choose j'. Equation (3.3) holds provided a , r 0 /2 (see §3.2 for a discussion). The first term in the numerator of equation (3.3) represents the effects of the viscosity mismatch, the second term embodies the effects of surface tension, while the final term shows the effects of activity.
As L n (0) ¼ 0 for all n, we obtain the following expression of s n (r 0 ) in the passive limit a ! 0: The first term (square brackets) of equation (3.5) is equivalent to the landmark result obtained by Paterson [58] for viscous fingering in a radial geometry (eqn. (10) in [58]), upon imposing the injection rate Q in Paterson's formula equal to the total amount of the droplet's growth per unit time in our setting (that is, Q ¼ pr 2 0 in dimensionless variables). The last term 1/(f þ 1), however, distinguishes our result from Paterson's, and stems from the fact that, here, the invading fluid also grows within the fingers. Figure 3 shows s n (r 0 ) versus n for various values of f and a, and for an unperturbed droplet radius r 0 ¼ 5. Only integer values of n (circles in figure 3) have a physical interpretation. The first mode n ¼ 1 corresponds to a translation of the droplet, and as L 1 (x) ¼ 0 for all x, s 1 (r 0 ) ¼ 0 for all values of r 0 , a and f. Higher modes n ! 2 correspond to the formation of n fingers on the interface of the active droplet and, if unstable (that is, if s n (r 0 ) . 0), could potentially initiate the multicellular protrusions observed in tumours [4].
In passive fluids, f . 1 (that is, the invaded fluid is more viscous than the invading one) is a necessary condition for fingering to be initiated, as instabilities can only grow when the pressure gradient near the interface is lower in the invading fluid [69]. However, we observe that s n (r 0 ) increases with a, so that modes that are stable when a ¼ 0 may become unstable in the presence of activity a . 0 (compare the dashed and solid blue curves, obtained with f ¼ 0.5, in figure 3). Therefore, activity can trigger fingering in systems that are stable otherwise, as well as enhance and/or change the dominant modes in droplets that are already unstable.
Activity lowers the pressure gradient of the invading fluid near the interface, hence promoting instabilities. Using the expression of p 0 and p 0 0 given in equations (2.5), we find that the condition for fingering, rpj r¼r0 , rp 0 j r¼r0 , is equivalent to f þ 2a/r 0 . 1 when the effects of surface tension are negligible.

High-activity regime
As already mentioned, jv s j ¼ a/b 10 29 m s 21 is a characteristic velocity at which the active fluid would move in an unbounded space, under uniform pressure and without growth. If jv s j exceeds the growthgenerated velocity at the interface, the active fluid's motion is frustrated and further instabilities occur across its entire area. We call this regime, for which dimensionless a . r 0 /2, 'high activity'. In this case, the derivation of equation (3.3) presented in the appendix, which assumes that perturbations are only arising at the interface and decaying away from it, is not valid. This regime is potentially relevant in the behaviour of real epithelial tissues, in which fingering at the boundaries is accompanied by velocity swirls forming across the entire area of the tissue [23].
Hence, equation (3.3) is only valid for the 'low activity' regime (a , r 0 /2). Yet, even in this regime, the active droplet may feature a region of instabilities near its centre r , 2a, where the velocity magnitude jv 0 j (given by equation (2.5a)) is less than jv s j. In particular, this situation would have also occurred in the history of the system considered in figure 3. In practice, a separate simulation-based study would be most appropriate to obtain the velocity field throughout the whole active region and to further examine the high-activity regime. In further sections of this paper, we only examine the system's behaviour in the low-activity regime.

Onset of fingering
For small enough radii, surface tension stabilizes the active droplet, but its strength decreases as that droplet grows. Therefore, there exists a critical radius for the onset of fingering, below which the active droplet grows circular and unperturbed. As n ¼ 2 is always the first mode to become unstable, r c (a, f) defined by the conditions s 2 (r 0 ¼ r c ) ¼ 0 and @ r0 s 2 (r 0 )j r0¼rc . 0 provides an estimate of that critical radius. Using equation (3.3) with L 2 (x) ¼ x(3x 2 4)/(2x 2 3), these conditions are equivalent to finding a polynomial root (first condition) within a subdomain (second condition), and have a unique positive solution.
We plot r c (a, f ) versus a in figure 4, which shows that increasing activity decreases the minimum radius for fingering. When f , 1 (blue curve in figure 4), there exists a minimum value of a, below which fingering cannot occur, because the higher viscosity of the invading droplet has a stabilizing effect. When f ¼ 1 (green curve in figure 4), a moderate increase of a may decrease r c multiple times. The impact of a in fingering is, however, reduced when f . 1 (the red curve in figure 4), because in that case the interface would be unstable even without activity.
The range of dimensionless r c , presented in figure 4, would correspond to a radius of 10-100 mm; however, it can be much higher for lower values of a. This range of r c is nevertheless in qualitative agreement with the tumour size at which the onset of fingering occurred in experimentally studied carcinomas [4].

Dominant mode
We now address the question of how many fingers are visible in practice, or, technically speaking, the question of which perturbation mode dominates during growth. Viscous fingering studies suggest  the dominant mode is the one satisfying the so-called maximum-amplitude criterion [70]. The criterion is satisfied by the mode n d experiencing the largest total aggregated growth in amplitude z n over the entire history of the system [70]. Following [70], we obtain z n by integrating the rate of the perturbation's growth, as predicted by our linear stability analysis, over that history, where R n is the radius at which mode n is first destabilized (i.e. the minimum radius at which s n (R n ) becomes positive). The dominant mode n d is then obtained for each r 0 from the conditions, @ n z n (r 0 )j n¼n d ¼ 0 ( 3 :7a) and @ 2 nn z n (r 0 )j n¼n d , 0, (3:7b) used to locate the maximum aggregated growth. As we shall see, the selection of the dominant mode depends on the particular kinetics of the tumour growth. Some experimental studies have shown that an initially exponential growth [71] (corresponding to a constant k) subsequently slows down with time (implying a decrease in average k) as the tumour enlarges and its resource supply becomes a limiting factor [72][73][74]. Other kinetics have been measured for various tumours and phases of growth, including sigmoidal regimes in which the growth stalls [72 -74]. The kinetics where the tumour's radius grows linearly with time also naturally emerges when the tumour proliferates only within an outer rim [74]. Such growth can also occur as a temporary feature in a sigmoidal kinetics.
We thus proceed to discuss in detail the selection of dominant modes in two kinetic models of tumour growth: an exponentially growing tumour, where k is uniform and independent of time and r 0 (t) ¼ r i e k(tÀti)=2 (from equation (2.1b) at the interface, and with dimensional variables; r i being the initial radius at time t i ); and a tumour with a radius growing linearly with time, r 0 (t) ¼ r i þ n(t 2 t i ), with n the constant and uniform velocity of the unperturbed interface. In the latter case, the growth rate k appearing in equation (2.1b) evolves with time.

Exponential growth
The integral given in equation (3.6) cannot be expressed analytically for all values of f and a, and we plot in figure 5 the relationship n d (r 0 ) obtained from numerical evaluation. When f . 1, and for low radii close to the onset of fingering, we observe that the activity has only a moderate influence on the selection of the dominant mode. For later growth, when r 0 becomes large, the viscosity mismatch is the governing cause of fingering and the activity plays no role in the selection of the dominant mode. We numerically observe the power-law variations n d / r 3=2 0 , independent of a. When f ¼ 1, higher activities promote the selection of higher modes. We obtain numerically, and for large r 0 , the scaling n d / r power law is independent of a. When f , 1 and a is sufficiently high, some low-n modes will become destabilized. However, these perturbations will re-stabilize and decay as r 0 increases further, because the stabilization from viscosity mismatch dominates as the radius of the droplet grows: active terms of equation (3.3) vanish when r 0 ! 1, while terms involving f remain constant in this limit. We may recover analytically the observed scalings for large n d and r 0 , and for f ! 1. We derive in appendix B.1 the following results when r 0 ! 1: n d % (f À 1) 1=2 Â r 3=2 0 for f . 1 and n d % [(a=2) 2 þ (a=2) 1=2 ] 1=2 Â r 1=2 0 for f ¼ 1, which we give below in dimensional variables to highlight the influence of the various physical parameters: when r 0 ) ð 2g bk Þ 1=3 , and with c % 0.06 defined as the smaller of the two solutions to 3c ¼ 3 þ lnc.

Linear growth
We also investigate pattern selection for a tumour with a radius growing linearly with time. The derivation of s n (r 0 ) proceeds along identical lines, although in this case, the integral in equation (3.6) can be expressed analytically. Details of this calculation are given in appendix B.2, and in this case we find that n d % c(a þ f 2 1) 1/2 Â r 0 is valid for all values of the physical parameters in the low activity regime when r 0 ! 1, and where c % 0.06 is the constant defined previously. We thus write, in dimensional form, when r 0 ) ð g bn Þ 1=2 (note that the dimensionless variables are defined differently in the linear growth, as explained in appendix B.2). The difference in the n d versus r 0 power-law dependency between equations (3.8) and (3.9) highlights the role of the growth kinetics in the fingering pattern, and is discussed with more details in the following.

Comparison and discussion
We now examine the evolution of the tumour's shape in the two growth kinetics studied above. An initially circular droplet is allowed to evolve, with the n-mode perturbation starting when the radius r 0 reaches R n , and with an initial amplitude of 0.2 (corresponding to approx. 2 mm). The perturbation is subsequently allowed to grow according to equation (3.6), such that z n (r 0 ) represents the weight of the n-mode at the unperturbed droplet radius r 0 . We further assigned a random phase difference between each n-mode perturbation.
We present in figure 6 examples of droplet patterns obtained with this procedure, where fingering is driven by either viscosity mismatch (left) or by activity (right), in both the exponential (top) and linear (bottom) growth regimes. In the linear growth, activity-driven fingers emerge more distinctively than in the passive droplet; the opposite is observed in the exponential growth. These results, as well as the analytical scalings presented above, demonstrate that the role of activity in fingering depends on the kinetics of the tumour's growth, and is indeed enhanced in the slower, linear growth kinetics. This assessment could potentially provide a basis for the mechanism behind the onset of metastasis, when the bulk growth of the primary tumour slows down or saturates.
The results presented in figures 5 and 6 relate direct observables of the tumour's geometry, and such measurements should indeed be envisaged by experimentalists. Note, however, from figure 6 that the number of fingers (approx. 10) visible at r 0 100 mm when f ¼ 1 is in agreement with the experimental observations shown by Cheung et al. [4].

Conclusion
We have devised a model of a growing and self-propelled tissue that isolates the role of four mechanical parameters (summarized in table 1) on its dynamics. The theory is based on experimental evidence and is analytically trackable. We used it to describe the evolution of an embedded two-dimensional circular rsos.royalsocietypublishing.org R. Soc. open sci. 5: 181579 droplet that could model a carcinoma in an epithelial layer. In this example, we were able to highlight the basic mechanical conditioning required to form interfacial instabilities, reminiscent of the classical viscous fingering, and that could explain the tumour protrusions observed at the onset of metastasis. We notably find that the tissue's active traction and growth kinetics are central to shape the instabilities' pattern and evolution.
Our model, and the example of its application presented here, could further help predict the minimum tumour size for metastasis, as well as the number of subsequent invasive fingers emerging from the initial mass. To the best of our knowledge, these observable geometric quantities have yet to be measured systematically in experimental studies.
The relative analytical simplicity of our model allows the investigation of more complex settings, such as heterogeneous tumours where active forces and/or growth are not uniform, or processes where these parameters are evolving with time or are dependent on one another. It also offers constitutive equations that can be used in simulations, and we envisage such studies for systems with high traction forces, where active motions are faster than the growth velocity, and which may indeed be relevant in aggressive forms of cancer. Acknowledgements. The authors thank Drs J. Prost, S. Lira, A. Hallou and P. Szymczak for their insightful comments.

Appendix A. Linear stability analysis
We here follow the lines of the demonstration given by Paterson [58]. We substitute the linearized perturbations to the hydrodynamic fields in both fluids, (dp, dp 0 , dv, dv 0 ), into equations (2.1) -(2.2) and obtain, to linear order @ r dp ¼ Àdv r , ( A 1 a) v r @ u dp ¼ ardv u À rv r dv u (A 1b) in the active fluid, and @ r dp 0 ¼ Àfdv 0 r , ( A 2 a) @ u dp ¼ Àrfdv 0 in the passive fluid.
Introducing the ansatz of periodicity in u, that is (dp, dp 0 , dv, dv 0 ) / e inu , in equations (A 1) -(A 2) allows us to calculate the partial derivatives with respect to u. Combining the resulting equations, and using the expression of v 0 of equation (2.5a), leads to second-order differential equations for the pressure perturbations in both fluids, n 2 dp ¼ (r À 2a)@ r (r@ r dp) ( A 3 a) and n 2 dp 0 ¼ r@ r (r@ r dp 0 ): (A 3b) We now assume a , r 0 /2 (see §3.2). Equations (A 3) are both solved by linear combinations of two functions. One function in each of these combinations has incorrect asymptotic behaviour (diverging, instead of decaying away from the interface) and is removed. The retained functions serve to formulate the allowed forms of pressure perturbations, dp ¼ q(t) X n j¼0 (À1) j n þ j n þ j j n j The functions q(t) and q 0 (t) are then obtained by finding dv r (using equation (A 4a) in equation (A 1a)) and dv 0 r (using equation (A 4b) in equation (A 2a)), and substituting the resulting expressions into the kinematic boundary conditions, equation (2.3). One obtains the following: dp ¼ r 0 n @ t dr À dr 2 1 þ 2a r 0 À L n 2a r À (n À 1) 2a r (A 5a) and dp 0 ¼ fr 0 n @ t dr þ dr 2 r r 0 with L n the function defined by equation (3.4), and where the u-dependency e inu has been incorporated into dr. Using both expressions in the pressure boundary condition, equation (2.4) linearized to first order in dr, and substituting dr / f n (t) e inu leads to an expression for @ t f n /f n . Upon using the definition of s n (r 0 ), equation (3.2), we finally obtain equation (3.3).
with A n ¼ fÀ1 n(nþ1) r 3 0 and B n ¼ ðnÀ1ÞðfÀ1Þ 3(fþ1) , by using the expression of R n obtained from solving s n (R n )j a¼0 ¼ 0. The dominant mode n d is calculated for each r 0 from the conditions of equations (3.7). Taking the limit r 0 ! 1, and reintroducing dimensional variables, lead to equation (3.8a).
When f ¼ 1, the selection of the dominant mode depends on the activity a. We first Taylor expand L n (x) to first order in x, L n (x) % 2n(n À 1) 2n À 1 x: This term grows linearly with n for n ! 1, while higher order coefficients in x plateau in this limit. Therefore, this expression of L n is efficient even for x 1, and we use it to approximate R n % ðnþ1Þð2nÀ1Þ 4a h i 1=2 and s n (r 0 ) % n(nÀ1)

B.2. Linear growth
In this kinetics, the growth rate k is not a constant, as opposed to the velocity n of the unperturbed interface which is now the adequate growth parameter. We thus use a different set of dimensionless variables, based on n. The characteristic length is now ' ¼ ( g bn ) 1=2 , and times are re-scaled by '/n, activities by a Ã ¼ bn, pressures by p Ã ¼ b'n and velocities by n. The derivation presented in appendix A proceeds similarly, and we find that the rate of perturbation's growth s n (r 0 ) in this system of variables is written as s n (r 0 ) ¼ 1 r 0 (f À 1)(n À 1) À n(n 2 À 1)=r 2 0 þ L n (a) f þ 1 þ L n (a) À na : The integral of equation (3.6), with the change of variables now written dt ¼ dr, may then be calculated analytically and we obtain the following: