The complex ecosystem in non small cell lung cancer invasion

Many tumors are characterized by genetic instability, producing an assortment of genetic variants of tumor cells called subclones. These tumors and their surrounding environments form complex multi-cellular ecosystems, where subclones compete for resources and cooperate to perform multiple tasks, including cancer invasion. Our recent empirical studies revealed existence of such distinct phenotypes of cancer cells, leaders and followers, in lung cancer. These two cellular subclones exchange a complex array of extracellular signals demonstrating a symbiotic relationship at the cellular level. Here, we develop a computational model of the microenvironment of the lung cancer ecosystem to explore how the interactions between subclones can advance or inhibit invasion. We found that, due to the complexity of the ecosystem, invasion may have very different dynamics characterized by the different levels of aggressiveness. By altering the signaling environment, we could alter the ecological relationship between the cell types and the overall ecosystem development. Competition between leader and follower cell populations (defined by the limited amount of resources), positive feedback within the leader cell population (controlled by the focal adhesion kinase and fibronectin signaling), and impact of the follower cells to the leaders (represented by yet undetermined proliferation signal) all had major effects on the outcome of the collective dynamics. Specifically, our analysis revealed a class of tumors (defined by the strengths of fibronectin signaling and competition) that are particularly sensitive to manipulations of the signaling environment. These tumors can undergo irreversible changes to the tumor ecosystem that outlast the manipulations of feedbacks and have a profound impact on invasive potential. Our study predicts a complex division of labor between cancer cell subclones and suggests new treatment strategies targeting signaling within the tumor ecosystem.


Introduction
Lung cancer is the second most prevalent type of cancer causing over 150,000 deaths per year in the United States [1]. Insufficient progress has been made in achieving efficacious treatments. One of the main barriers in developing new treatment strategies is the vast diversity between and within cancers; heterogeneity exists between patients with the same tumor type, between tumor loci within a patient (i.e. metastases and primary tumor), and within the primary tumor itself [2,3]. Cancer is distinguished by loss of normal control over cell processes leading to genetic instability and unregulated growth. Genetic instability creates array of different clonal populations with different cell fitnesses, renewal and invasion potential [4]. Competition between different cancerous subclones and between cancerous and normal cell types sets the stage for classical ecological dynamics in the tumor microenvironment. The outcome of this process determines success of the tumor progression and its understanding may help discover novel treatment strategies [5,6].
Invasion of surrounding tissue, either locally or distally via metastasis, is a hallmark of cancer [7]. Extensive research has detailed that invasion is mediated by interactions between tumor and extracellular matrix [8,9] and cancer-associated fibroblasts [10], but there is a lack of focus on the cooperative interactions between different cancer cell types, either phenotypically or genotypically distinct. Indeed, in mouse models of lung cancer, collective invasion of cancer cells was shown to correspond markedly more successful metastasis [3,[11][12][13], confirming the critical role of collective invasion in driving cancer progression.
We recently developed a novel image-guided genomics approach termed SaGA that allowed us to identify at least two distinct phenotypic cell types in lung cancer invasion packs: highly migratory leader cells and highly proliferative follower cells [14]. Genomic and molecular interrogation of purified leader and follower cultures revealed differential gene expression prompting distinguishing phenotypes. Specifically, leader cells utilized focal adhesion kinase signaling to stimulate fibronectin remodeling and invasion. Leader cells also overexpressed many components of the vascular endothelial growth factor (VEGF) pathway facilitating recruitment of follower cells but not the leader cell motility itself [14]. However, leader cells proliferated approximately 70% slower than follower cells due to a variety of mitotic and doubling rate deficiencies. These deficiencies could be corrected by addition of cell media extracted from the follower only cell cultures, leading to conclusion that follower cells produce an unknown extracellular factor responsible for correcting mitotic deficiencies in the leader cells. In sum, leader cells provide an escape mechanism for followers, while follower cells (and follower cell media only) help leaders with increased growth. Together, these data support a service-resource mutualism during collective invasion, where at least two phenotypically distinct cell types cooperate to promote their escape.
In this new study, we developed population-level computational model to explore impact of the complex interactions between leaders and followers cell types on cancer progression. The model implemented effects of critical signaling factors controlling the communication between cell types and the interaction between cells and environment. We derived analytic boundaries dividing parameter space, representing the major signaling feedbacks, by the critical changes to invasion dynamics. Our study predicts the critical role of specific signaling pathways involved in the symbiotic interactions between cancer subclones for the overall success of cancer progression.

Methods
Our model tracks the cell counts of leader cells, L, and follower cells, F, the concentrations of extracellular factors VEGF, V, an unidentified Proliferation factor, P, and Fibronectin, N, as well as the size of the domains for leader cells, O L , and for follower cells O F . Based on the available data [14], the following processes have been implemented. Leader cells can expand their domain, O L , by secreting Fibronectin, which in turn relaxes competitive pressure on leader cell growth. Leaders also secrete VEGF, which is taken up by follower cells and causes follower cells to follow them. This was modeled by increasing the domain for follower cells, which in turn relaxes competitive pressure on follower cells. Follower cells secrete an unknown proliferation signal that increases the reproductive capacity of leader cells (initially smaller than follower cells). Leader and follower cells also must compete with each other for resources at rate c (see Fig 1).
We modeled cell counts (L and F species) as standard Lotka-Volterra competition [15]. The carrying capacity of the leader cells was dynamic and dependent on the amount of proliferative signal, P, present. This capacity increased in a saturating manner with P, with maximum equal to the follower cell carrying capacity, K F0 . Intra-and inter-specific competition was driven by Here r L and r F denote the rate of expansion for leaders and followers, respectively. The parameter c denotes the strength of competition between the two cell types-the degree to which the presence of one will impact the capacity of the other [16,17]. The capacity of the environment for follower cells is given by the parameter K F . The capacity for leaders depended on an initial capacity, K L0 , and on the amount of proliferation signal in a Hill-like manner with EC50, δ. Each extra-cellular species (V,P,N) had a production rate, β, and a degradation rate γ, the domain size variables (O L and O F ) also had a parameter denoting initial capacity (O L0 and O F0 ).

Reduction and feedbacks
Previous 3D spheroid experiments show that invasion occurs on a much faster time scale than reproduction [14]. By assuming that factors (V,P,N) and domains (O L ,O F ) change much faster than cell counts (equivalently γ V ,γ P ,γ N ,γ OL ,γ OF ) r L ,r F ), one can reduce these equations to a set of two equations (L,F), where variables in Eqs (3)-(7) are at their equilibria Using this reduction drastically decreased the complexity of the system. This allowed a twodimensional phase-space representation, facilitating the presentation of the results, and analytical derivation of many of the bifurcation conditions in the system.
To perform this reduction, we first defined the feedbacks based on the reduced system. The feedback that determines the leaders impact on their own domain expansion was denoted by , for the strength of the leader only feedback. The feedback that determines the leaders impact on follower cell growth was denoted by s LF for the strength of the leader to follower feedback. The feedback that determines the followers impact on leader cell growth was denoted by s FL ¼ b P g P d , for the strength of the follower to leader feedback. Second, using these assumptions, we re-wrote the leader-follower system as where Using this reduction we can derive several critical points in invasion. The reduced system  Table 1.

Transcritical bifurcation at zero
To determine the critical points in the leader-follower system, we calculated the Jacobian of the reduced system evaluated for the leader extinction equilibrium (O 1 : Here 1þs FL ÁF LE , the value of K L when F = F LE . The Jacobian has

Parameter(s) Value(s) Interpretation
Variable Feedbacks The strength of the feedback within leader cell population The strength of the feedback from leaders to followers For c <

Saddle node bifurcation
The system undergoes a saddle node bifurcation when two coexistence equilibria (O 4 and O 5 ), representing non-zero populations of both leaders and followers, coincide and disapper. Beyond this bifurcation point the leader/follower populations undergo unbounded growth. This bifurcation was determined numerically using MatCont [18]. We found that this bifurcation point depends critically on both the leader feedback strength, s L , and on the competition strength, c. One of these coexistence points is effected by the transcritical bifurcation, below.

Transcritical bifurcation at infinity
When the leader feedback strength is sufficiently high relative to competition, leaders and followers may undergo unbounded growth from the initial conditions belonging to the certain regions of the phase space. We describe this scenario as an attractor basin in the phase space for the stable infinity attractor. However, if s L is reduced (or c is increased) beyond a certain threshold, infinity becomes unstable. This corresponds precisely with the loss of an unstable coexistence equilibrium with non-zero values of both leaders and followers (O 5 ). Leaders and followers that are coexisting must satisfy In the case that follower populations are large relative to δ, which has a discontinuity at defining the loss of one of the coexistence equilibrium points (O 5 ) when it moves to infinity. We describe this as the transcritical bifurcation at infinity as the stability of infinity changes at this point.

Leader and follower ecosystem
Leader and follower cell types in non-small cell lung cancer spheroids were previously isolated using a fluorescence technique termed SaGA [14] ( Fig 1A). We found that leaders and followers are genotypically and phenotypically distinct populations of cancer cells that exchange a variety of signaling molecules to coordinate complex behavior during invasion. In this new work, we focus on four main channels of communication (see Fig 1B). Through the activation of focal adhesion kinase (FAK), leader cells secrete fibronectin in an autocrine manner. This leads to ECM restructuring and expansion of leader cell domain, O L , (see Methods) which ultimately increases the leader cell count. The strength of this positive feedback is characterized in our model by s L (strength of Leader only feedback). Leader cells also secrete VEGF. In the leader-follower ecosystem this promotes follower cells to track expanding leader cells, increases follower domain size (O F ), and ultimately, follower cell count. In our model, the strength of this feedback is given by s LF (strength of Leader to Follower feedback). Follower cells secrete an undetermined proliferation signal, as evidenced by the observation that follower-only cell media increases leader cell growth rate [14]. The strength of this feedback is given by s FL (strength of Follower to Leader feedback) in the model. Finally, both cell types compete for the same resources. This has an effect of limiting the capacity of the each cell type through competition, modeled here by the feedback c [16,17]. These feedback mechanisms were incorporated into a modified Lotka-Volterra type competition-cooperation model. We chose a Lotka-Volterra model to focus on the ecological aspects of competition in the cancer ecosystem. Here, the leader cells could grow to a total capacity K L , which is an increasing function of the proliferation signal secreted by the follower cells. This capacity was reached when a combination of leader and follower cell densities (cell counts divided by domains) exceeds K L (see Methods). Increases in the domain size of each type (by Fibronectin secretion in the leader case and VEGF in the follower case) limited the overall density of that cell type and mitigated its impact on the overall capacity of the system. Increasing competition, for example by limiting resources, increased the impact of either cell type on the conjugate capacity type (e.g. how leader density, L/O L , impacts follower capacity K F ).
This system of the feedbacks between the leader and follower cells describes a complex dynamical ecosystem. The impact these feedbacks may have on cancer growth or invasion is unclear. Leader and follower cells are engaged in competition for resources but can also be engaged in cooperation and play supportive roles. For example, invasive leader cells provide new territory for the follower cell population and are supported by proliferative follower cells. In the following, we analyzed the model to find critical turning points for the ecosystem dynamics.

Multiple types of invasion dynamics
We found that multiple feedbacks between the leader and follower cell populations could produce a wide variety of complex dynamics. When competition strength, c, was high and the strength of the leader only feedback, s L , was moderate, population dynamic was bounded and resulted in a stable cell count for both leader and follower cell populations, with the former decaying to zero, as well as a stable domain size (Fig 2A). In contrast, when feedback was large and competition was moderate, population dynamics revealed an unbounded growth ( Fig 2B). Intermediate values of both c and s L led to dynamic regimes that depended on the initial cell count: ecosystems with large initial cell count underwent unbounded growth, while small ecosystems attained a stable size ( Fig 2C). Many of the predictions of our study are based on analysis of a reduced model system we assumed that extra-cellular factors, such as VEGF and fibronectin, change much faster than leader and follower cell counts. We found that dynamics of leaders and followers were consistent with the full system across many orders of magnitude of the ratio of these two rates (Fig 2D). These types of dynamics are in a qualitative agreement with experimental studies which revealed (a) rapid expansion of intact leader-follower ecosystem and (b) that blocking specific feedback mechanisms in vitro can reduce or block cell population growth. Specifically, blockade of fibronectin signaling or blockade of VEGF signaling led to significantly reduced invasion [14].
This array of behaviors can be explained by the critical shifts in the cell population dynamics due to the changes in the feedbacks strength. We found that depending on the level of competition, c, and the strength of invasiveness of leaders, s L , the leader-follower ecosystem can operate in one of five different regimes, as described below (Fig 3).
Leader extinction. When competition was high and invasive feedback was minimal, the leader cells (the weaker competitor) were forced to extinction while the follower cells persisted and its population reached a stable size (Fig 3A and 3B). There was a critical level of competition between leaders and followers, given by c > K SS L K F (see Methods, Transcritical bifurcation at zero for derivation), required for this type of dynamics. This critical level of competition, the ratio of the capacity of leader cells, K SS L , to that of the follower cells, K F , essentially depends on the fitness differences between leader only and follower only cell populations. Leader and follower populations with similar fitness would tolerate a much higher competition threshold without driving one species to extinction. From the dynamical systems perspective, when c > K SS L K F and s L is sufficiently low (see below), the only stable equilibrium in the phase space is O 1 and all system trajectories converge to this equilibrium point representing the leader cells extinction state (Fig 3B).
Leader extinction with escape. If competition was above the leader extinction limit, c > K SS L K F0 , but not high enough to balance the impact of the leader only feedback, c < K F s L − 1, there were two possible outcomes depending on the initial population size (see Methods, Transcritical Bifurcation at Infinity for derivation) (Fig 3A and 3C). The second condition, c < K F s L We define ε = r F /γ and show simulations of the full system with ε between 0.001 and 1 for leader cell dynamics with parameter values as in C. There was less than a 5% difference in growth and decay rates for ε < 0.1. https://doi.org/10.1371/journal.pcbi.1006131.g002 The complex ecosystem in non small cell lung cancer invasion − 1, can be interpreted as a balance between positive feedback, s L , and negative feedback, c. In this regime, leaders could go extinct if the initial population of leader cells was sufficiently small. Alternatively, if initial populations of leaders and followers both were large enough, the ecosystem could grow unboundedly. Thus, our model predicts, that the ability to undergo successful collective invasion depends on whether the initial bulk size is larger than a critical amount. These types of dynamics with divergent outcomes occur when competition is large enough to be able to drive leaders extinct, but small enough so it can be outbalanced by the strong invasive effects of the leader cells.
In the phase space of the model, the basins of attraction of the two distinct dynamical regimes are separated by a critical boundary (separatrix of a saddle equilibrium O 5 ) where the cell bulk size determines its ultimate fate (Fig 3C). Both infinity and the leader extinction equilibrium (O 1 ) are stable attractors representing two possible end solutions of the system dynamics.
Non-invasive dynamics. When competition was smaller than the extinction limit, c < K SS L K F0 , but large enough to balance leader feedback strength, c > K F s L − 1, the ecosystem size remained bounded and both leaders and followers attained a stable and non-zero population size. In the phase space, this type of dynamics corresponds to conversion to the stable equilibrium O 4 (Fig 3A and 3D). We refer to this as non-invasive dynamics, as the cells cannot grow beyond a defined size. In this case, while competition was present, it was too weak to lead to extinction, while leader population was not invasive enough to promote unlimited growth. This scenario represents stable, non-invasive dynamics.
Multimodal dynamics. If competition was (a) small enough to allow leader existence, c < K SS L K F0 , (b) small enough relative to the leader feedback strength, so that escape was possible, c The complex ecosystem in non small cell lung cancer invasion < K F s L − 1, but (c) high enough, so that for small initial population of leader it could balance the positive leader feedback, leader and follower cell dynamics depended on the initial population size (Fig 3A and 3E). Ecosystems with a large initial cell count would grow without bound but those with a small initial cell count would reach a stable population size, due to the competition as in the non-invasive dynamics case. On the phase plane the last outcome was represented by a contraction to a stable equilibrium O 4 . This critical boundary was defined by a separatrix of a saddle fixed point O 5 (Fig 3E). (This separatrix was determined numerically by reversing time [19].) Aggressive dynamics. When leader invasive strength was sufficiently high and competition was sufficiently low, the only possible outcome was unbounded growth of both cell populations (Fig 3A and 3F). In this case, the only stable attractor in the phase space is infinity where all system trajectories are converged to.
In summary, our analysis revealed that the complex balance of the feedbacks in the leaderfollower ecosystem can lead to the multiple types of population dynamics. When the leaders' invasiveness was low, the outcome depended on the competition between two populationsstrong enough competition promoted leader extinction, while weak competition allowed stable coexistence states with bounded size of both leader and follower cell populations. As leader invasiveness rate increased, the system revealed a new state with unbounded growth. This aggressive dynamic state coexisted with a stable attractor representing a bounded size of both populations if competition between leaders and followers was strong enough. Otherwise, unlimited population growth was the only outcome. Based on the system dynamics derived above, next we will show how critical boundaries between parameter regimes could be exploited to lead to profound changes in the ecosystem dynamics.

Limiting leader feedback leads to irreversible changes in invasion
During multimodal dynamics (e.g., Fig 3E), leader and follower cell populations can undergo explosive growth or achieve a stable count depending on the initial size of the ecosystem. We examined the impact of limiting the invasive leader feedback in scenarios of this type (Fig 4). Even when the ecosystem was initially sufficiently large to support unbounded growth, after reducing invasive leader feedback s L (Fig 4A), the ecosystem was forced into the non-invasive dynamics type and the total bulk of the cell population reduced reaching a steady-state ( Fig  4E). Importantly, the leader and follower cell populations remained stable and bounded after restoring invasive leader feedback to its original strength (Fig 4E, right side).
From the point of view of the dynamical systems analysis, reducing leader feedback changed the phase space, so the only stable attractor was non-zero equilibrium (O 4 ) (Fig 4C). In this regime, unlimited growth was abandoned and the system converged to the equilibrium state (O 4 ) corresponding to the bounded size of both cell populations. This equilibrium remained stable even after the feedback was restored to its original level (Fig 4D).
Our model predictions (Fig 5A) are consistent with in vitro data (Fig 5B). Using siRNA blocking we previously showed that expression of fibronectin (which is characterized by the strength of leader only feedback, s L , in the model) led to the low invasion potential and a stable cell population size [14].

Increasing competitive signals leads to leader extinction
We next tested effect of increasing competition between leader and follower cell populations on the ecosystem dynamics (Fig 6). Leader cells excrete extracellular factors that induce the death of the followers and leaders alike [14], which supports competition. Here, we again started from aggressive unbounded type of dynamics and then increased competition strength ( Fig 6A). This caused change of the ecosystem dynamics. Both cell populations reduced the size, with leader cell population going to extinction state ( Fig 6E). However, upon restoring   [14]. C) In the model, blocking leader to follower feedback, s LF , limited invasive area and cell count. D) Impact on invasion of leader cell cultures during siRNA block of VEGFR2 (siKDR). Scale bar 100μm. Invasive area was significantly reduced after VEGFR2 block (p<0.0001) (right), compare to control (left). Reproduced from [14]. https://doi.org/10.1371/journal.pcbi.1006131.g005 The complex ecosystem in non small cell lung cancer invasion competition to the original level, leader and follower cells reemerge and grow unboundedly again. The last can be avoided if no leader cells remain (complete extinction).
Again, this dynamic can be easily understood using bifurcation analysis. Increasing competition strength made leader extinction equilibrium state O 1 stable (Fig 6C). However, when competition was restored to its original level, O 1 became unstable again and leader and follower cells returned to escape dynamics ( Fig 6D). Importantly, in the extreme case of very small cell populations, cells undergo discrete and stochastic dynamics and complete extinction of a small population of leaders is possible in a finite time, leading to irreversible changes due to competition increase (similar dynamics was described in our previous study [20]).

Support for leaders has large impact on aggressiveness
Changing the strength of the feedbacks that determine the interaction between leaders and followers (s LF and s FL ) could also impact the dynamics. Leader cells secrete VEGF (denoted here by s LF ) that helps follower cells to expand their territory and follower cells secrete a proliferation signal (denoted here by s FL ) that allows leaders to increase their proliferative capacity. These two feedbacks have distinct impacts on the overall ecosystem dynamics. Perturbations to s LF (changing the impact that leaders have on followers) changed the system dynamics (assuming that the cell count was small enough at the time of the intervention) from unlimited growth to the bounded type. The size of both leader and follower cell populations decreased reaching non-zero steady-state (Fig 7E). This regime persisted as long as the feedback from the leaders to followers remained low. However, increasing s LF to its original level restored the system dynamics with unlimited cell population growth (Fig 7E, right side).
Using bifurcation analysis, we found that reducing impact that leaders have on followers shifted the location of the saddle node bifurcation boundary that separated state with The complex ecosystem in non small cell lung cancer invasion unlimited growth only dynamics and a state with coexistence of the unlimited growth and a stable equilibrium attractor (O 4 ) regimes ( Fig 7A). Effectively, decreasing s LF increased the threshold level of the invasive leader feedback (s L ) needed to cause unbounded growth. Thus, reducing s LF made the system to converge to the stable equilibrium state O 4 corresponding to the bounded size of both cell populations (Fig 7C). However, increasing s LF to its original level changed the phase space again, so infinity became the only stable attractor (Fig 7D) and unlimited growth dynamics resumed.
Our model predictions (Fig 5C) are consistent with in vitro data ( Fig 5D). Using siRNA to block the VEGF receptor VEGFR2 (siKDR in Fig 5D), we previously showed that blocking the leader to follower feedback led to the limited invasion potential and stable cell population size (Fig 5D) [14].
Finally, we tested the role of the follower to leader feedback (s FL ) and found that perturbations to s FL have a significant impact on the system dynamics. In contrast to s LF , changes to the s FL changed both the location of the saddle node bifurcation boundary and the transcritical bifurcation boundary of the leader extinction ( Fig 8A). Therefore, decreasing s FL both increased the threshold on the leader invasion strength (s L ) needed to cause unbounded population growth and decreased the threshold of the competition strength (c) needed to induce leader population extinction. We have exploited this to show that decreasing s FL can cause irreversible change in the cell population bulk. Again, starting with unlimited growth dynamics (Fig 8B), decreasing follower to leader feedback, s LF , reversed the dynamics and both leader and follower cell population reduced in size converging to the steady-state (Fig 8E). This regime with bounded ecosystem size persisted after the feedback was restored (Fig 8E, right  side). The complex ecosystem in non small cell lung cancer invasion Using dynamical systems analysis, we found that reducing follower to leader feedback (s FL ) triggered the system convergence to the stable attractor (O 1 ) representing the leader extinction state (Fig 8C). When the feedback was restored, O 1 becomes unstable but the ecosystem fell to the attraction basin of the stable equilibrium O 4 and avoided regime of unlimited growth ( Fig  8D). In a more general case, the outcome depended on the balance between the leader to follower, s FL , and follower to leader, s LF , feedbacks, with higher s LF requiring more significant s FL decrease to avoid unbounded growth (Fig 8F).

Summary of perturbations to cancer ecosystem
A complex balance of the feedbacks within the cancer cell ecosystem allows for some alterations of the feedback parameters to have significant impacts on the ecosystem dynamics. We summarized these different possibilities in Table 2 from the perspective of achieving the goal to reduce cell population bulk. Hence, manipulating s L , s LF , s FL should be interpreted as decreasing these feedbacks, whereas manipulating c should be interpreted as increasing c. We also examined the possibility of non-targeted cell death, such as might occur during non-specific chemotherapy (implemented via a non-targeted "enforced" reduction of the cell population). Manipulations were either irreversible, so the system dynamics remained altered upon cessation of the perturbation (e.g. irreversible leader extinction or irreversible stabilization of the cell count), or caused only temporal and reversible reduction of the cell bulk. In some cases, such as leader extinction with escape and multimodal dynamics (see Fig 3), the size of the initial cell bulk dictated possible outcomes of the feedback perturbations. The outcomes described in Table 2 represent the best-case scenario. Thus, perturbations were started from an appropriate initial state and maintained long enough to achieve the desired effect. This analysis revealed that certain parameter regimes are more sensitive to the perturbations than others. Specifically, in the leader extinction with escape regime (area (2) in Fig 3A) and the multimodal dynamics regime (area (4) in Fig 3A) perturbations could have irreversible impacts on the ecosystem. In these cases, any perturbation (death, reduction in s L , s LF , s FL , or increase in c) can potentially force the system to cross the critical boundary (separatrix) and transition from explosive growth to a steady-state dynamic. These regimes give a unique opportunity to impact the invasiveness of the ecosystem.
Also, certain perturbations could force the ecosystem into a state where leader extinction (O 1 ) is stable. This occurs when applying sufficient increases in the competition pressure, c, or decreases in the support from followers to leaders, s LF . In these cases, it is possible for the discrete and stochastic nature of the cell population dynamics to define the ecosystem fate. Thus, a sufficiently long perturbation could irreversibly eradicate a sufficiently small discrete number of leader cells [20].

Plasticity of leader and follower phenotypes
Leader and follower cells may have distinct phenotypes because of underlying genetic or epigenetic differences. For an example of the latter, during angiogenesis, stimulation of endothelial cells by VEGF creates the epigenetic emergence of tip and stalk cell phenotypes with distinct roles in new vessel formation [21]. It remains an open question as to which is the case with leader and follower cells in lung cancer. We have shown that: a) leader cells are a phenotypically distinct subpopulation of lung cancer cells; b) leader and follower cells show distinct genetic expression profiles [14].
Observations of leader only populations over multiple months show that the leader cell phenotype is stable over many generations, maintains invasive morphology, and does not return to the follower phenotype. Similar observations of follower only populations do show an emergence of leader cell phenotype suggesting that follower cells can convert to the leader phenotype. This would suggest that the emergence of leaders from followers is due to epigenetic plasticity and additional unpublished data from the Marcus lab support this.
In accordance with the possibility that follower cells can convert to the leader phenotype, we have completed simulations to show our main result-the relative strengths of invasiveness and competition determine the dynamic regime of the leader follower system-remains true in this case. We found that when follower cells can convert to leaders, leader extinction is no longer possible. However, in agreement with previous results we observed that: i) high

Leader Extinction w/ Escape
Infinity attractor basin c, s FL , s L , s LF , death Irreversible leader extinction.
Stable attractor basin c, s FL , s L , s LF , death Reversible cell bulk reduction.

Non-invasive Dynamics
N/A c, s FL , s L , s LF , death Reversible cell bulk reduction.

Multimodal
Infinity attractor basin c, s FL , s L , s LF , death Irreversible stabilization of cell bulk.
Stable attractor basin c, s FL , s L , s LF , death Reversible cell bulk reduction. The complex ecosystem in non small cell lung cancer invasion competition and low invasiveness led to only tumors of fixed size; ii) low competition and high invasiveness led to tumors growing unboundedly; iii) moderate competition and invasiveness led to scenarios where the outcome of tumor growth depended on initial conditions (Fig 9).

Discussion
Heterogeneity of tumors, at the genetic, epigenetic, and phenotypic levels, is one of the main obstacles to developing new effective treatment strategies. Tumor cells rapidly evolve forming highly efficient symbiotic systems with well-defined labor division targeted to augment tumor survival and expansion. In lung cancer collective invasion packs observed in vitro, two distinct populations of cancer cells-highly migratory leader cells and highly proliferative follower cells-have been recently identified [14]. In this new study, we used computational models to explore collective dynamics of the leader-follower ecosystem and to exploit approaches that can effectively disrupt it. We found that competition between two populations (defined by the limited amount of resources), the positive feedback within the leader cell population (controlled by the focal adhesion kinase and fibronectin signaling) and impact of the follower cells to the leaders (represented by yet undetermined proliferation signal) all had major effects on the outcome of the collective dynamics. While increase of the positive feedback within the leader cell population would ultimately lead to the system state with unbounded growth, manipulating follower to leader feedback or increasing competition between leader and follower cell populations was able to reverse this dynamic and to form a stable configuration of the leader and follower cell populations.
Our model highlights the importance of fibronectin remodeling in invasion. Fibronectin is the major ligand of the focal adhesion kinase (FAK) pathway. Our previous empirical work showed that FAK signaling was a key distinguishing feature between leader and follower cells and critical for invasive leader behavior [14]. While we do not model FAK directly, our model predicts that fibronectin remodeling is the main driver of invasion by leader cells and disruptions in the FAK driven feedback loop will cause critical changes in the leader-follower population dynamics. Indeed, FAK is a well-known regulator of the tumor micro-environment: promoting cell motility and invasion [22]. FAK expression is upregulated in ovarian [23] and breast cancer [24] tumors with expression levels correlating with survival [25,26]. Many FAK inhibitors, such as defactinib, are currently in clinical trials with promising results [22,27,[27][28][29][30][31]. A key advantage of FAK inhibitors is that they impact both the tumor itself and the surrounding stroma where tumor associated fibroblasts also utilize FAK signaling to promote tumor invasiveness [32,33].
While commonly associated with angiogenesis in healthy and cancerous tissue, our previous work showed that VEGF mediates communication between leader and follower cells [14]. There is a long history of targeting VEGF to limit tumor invasiveness [34,35]. While great success has been seen in preclinical models [36,37], only moderate success was seen in clinical trials with anti-VEGF drugs such as bevacizumab [38,39]. This is largely due to cancers developing resistance to specific VEGF-therapeutics. In our model, VEGF stimulated followers to shadow leaders and expand their domain. However, we found that inhibition of VEGF had little impact on the ecosystem dynamics relative to the perturbations of the other axes (such as FAK or competition for resources).
Competition for resources is one of the principal forces that structures any ecosystem, including tumor ecosystems [6,40]. Our modeling work predicts that competition was a critical component in the leader-follower ecosystem. We found that when the strength of competition exceeded a critical threshold, leaders (the weaker competitor) were driven to extinction. Further, enhancements of the competition in the model changed the fundamental cell population dynamics. In some cases this meant stopping unbounded growth and promoting the extinction of the leader cells. Our previous in vitro work demonstrated that leaders may inhibit the growth of followers through an unknown secreted factor in cell media [14]. While still in the early stages, exploiting this inhibition may also provide similar benefits to those shown here as increases in competition.
Our previous study also revealed a currently unknown extracellular factor secreted by followers that corrects mitotic deficiencies and enhances leader proliferation [14]. Our modeling highlights this factor as having critical impact on the ecosystem dynamics. We found that blockade of this proliferation factor, modeled here by the strength of the follower to leader feedback, can cause critical shifts in the population dynamics. More work needs to be done to identify and understand the mechanism of this action, but preliminary results suggests that this may be a potential novel treatment axis that specifically targets the mutualistic interaction between leaders and followers.
Ecological forces shape the exchange of biomaterial between different biotic and abiotic environmental agents. These forces determine capacity of the ecosystem for different species (subclones) and the environment ultimately sets the fitness of each of the competitors. Classic ecological theory dictates that an abundance of many similar species (such as similar subclonal populations) will lead to a high competition for resources [41,42]. This competition can force the exclusion of inferior competitors and ultimately may reduce heterogeneity of the system. However, when symbiotic and mutualistic interactions occur, otherwise competitive species support each other and increase the capacity of the ecosystem [43,44]. Symbiosis between different subclonal populations may be particularly important during critical times when the tumor survival is in peril (such as hypoxia, metastasis or therapy). One critical moment in tumor progression occurs when highly proliferative tumor cells saturate the resource potential of their current environment. In order to obtain more resources, tumors need to invade new territory.
Collective invasion is a spatial phenomenon, where leader cells form the invasive periphery trailed by follower cells to invade new territory. We used a simplified non-spatial model here to make possible a quantitative bifurcation analysis that allowed us to calculate critical shifts in cell dynamics in response to the shifts in the biophysical model parameters. We focused our efforts on the balance of symbiotic and competitive effects in this complex tumor ecosystem: largely non-spatial phenomenon. To capture the benefits of invading novel territory to the cancer system, we used dynamic variables to represent the domain size of both leader and follower cells (O L and O F ). The advantages of this approach include: (1) biophysical properties are lumped into a few effective parameters with clear biophysical meanings, e.g., strength of interaction between populations, providing qualitative understanding of the complex interactions; (2) the low dimensional parameter space allows for systematic analysis to explore and determine the critical boundaries corresponding to disruption of collective invasion. While these simplifications allowed for useful analytical techniques, they also come with some limitations.
Our model assumed that leader cells, follower cells, and extracellular factors (VEGF, Fibronectin, etc.) were distributed in a spatially homogeneous manner, which is likely not the case in vivo. In particular, for successful collective invasion, leader cells must be at the periphery of an invasive front. Further, VEGF acts as a chemo-attractant in healthy cells and in most cancers [21,34], stimulating other cells to move up the VEGF gradient. These concepts of invasive front or VEGF gradient cannot be captured in our non-spatial model. Our model also fails to describe the motion of cells that cannot be accounted for by a simple increase in domain size, e.g., the impact of leader cell invasion to free up more room for growth and to colonize new sites. Despite these limitations the model reproduces many essential properties of the complex interactions found in experiments with the leader-follower tumor ecosystem and makes predictions about critical tipping points in the collective invasion of simplified leader follower cell populations.
Future work is necessary to extend our results to models incorporating the spatial evolution of leader cells, follower cells, and extracellular factors (VEGF, Fibronectin, etc.). One possibility includes using cellular Potts models to study invasion in cancer [45] to derive a reaction diffusion simplification using the procedure outlined in [46]. This approach should produce a spatially dependent continuous probability density approximation of a discrete and stochastic model. This model could allow us to extend our understanding of the spatially important aspects of the tumor ecosystem dynamics addressed here as well as investigate novel phenomena such as the impact of the extra cellular matrix organization and the interactions with cancer associated fibroblasts on collective invasion. Previous results to model complex tumor cell population dynamics range from detailed cellular level models (e.g. [9,[47][48][49]) to continuous models with a different degree of complexity (e.g. [20,[50][51][52][53]) similar to that proposed in our new study. While cellular level models can directly incorporate heterogeneous cell types and intrinsic tumor properties, including proliferation, metabolism, migration, protease and basement membrane protein expression, and cell-cell adhesion, they typically have high-dimensional variables and parameter space that is difficult to explore. Advantages of the reduced type of models include the low dimensional parameter space, where parameters have clear biophysical meanings, and which allows for systematic analysis to rapidly explore and determine the sensitive parameter space. We previously applied this approach to study cell interactions in chronic cancers and predicted conditions for explosive tumor growth [20]. Similar approach was applied to model cancer cell population dynamics in many other types of cancer [50,53,54].
The vast diversity between different cancers and even between different cell types within a single tumor remains one of the biggest hurdles to overcome to achieve personalized cancer treatment. This diversity leads to a complex array of interactions between different tumor cell types and the healthy surrounding tissue: the tumor ecosystem. Our work has isolated phenotypically unique lung cancer cells and taken a dynamical approach to understanding the interactions within the tumor ecosystem. We identified the critical features and interactions composing the leader-follower ecosystem, to explore vulnerabilities of the lung cancer invasive cell populations.