Invasion front dynamics of interactive populations in environments with barriers

Invading populations normally comprise different subpopulations that interact while trying to overcome existing barriers against their way to occupy new areas. However, the majority of studies so far only consider single or multiple population invasion into areas where there is no resistance against the invasion. Here, we developed a model to study how cooperative/competitive populations invade in the presence of a physical barrier that should be degraded during the invasion. For one dimensional (1D) environment, we found that a Langevin equation as dX/dt=Vft+Dfη(t)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$dX/dt=V_ft+\sqrt{D_f}\eta (t)$$\end{document} describing invasion front position. We then obtained how Vf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_f$$\end{document} and Df\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_f$$\end{document} depend on population interactions and environmental barrier intensity. In two dimensional (2D) environment, for the average interface position movements we found a Langevin equation as dH/dt=VHt+DHη(t)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$dH/dt=V_Ht+\sqrt{D_H}\eta (t)$$\end{document}. Similar to the 1D case, we calculate how VH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_H$$\end{document} and DH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_H$$\end{document} respond to population interaction and environmental barrier intensity. Finally, the study of invasion front morphology through dynamic scaling analysis showed that growth exponent, β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}, depends on both population interaction and environmental barrier intensity. Saturated interface width, Wsat\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W_{sat}$$\end{document}, versus width of the 2D environment (L) also exhibits scaling behavior. Our findings show revealed that competition among subpopulations leads to more rough invasion fronts. Considering the wide range of shreds of evidence for clonal diversity in cancer cell populations, our findings suggest that interactions between such diverse populations can potentially participate in the irregularities of tumor border.

www.nature.com/scientificreports/ poor clinical outcomes. On the other hand, the notion that invasion front geometry might reveal the driving mechanism behind the invasion 34,35 has sparked various studies on scaling properties of cancer cells invasion front in vivo [36][37][38] and using different mathematical models [39][40][41][42][43] . Despite the development of a diverse range of models on tumor invasion, clonal interaction remained overlooked. Motivated by interactions for cancer cells and inspired by a model on cooperative populations in the presence of environmental barriers 44 , we developed a model to study how environmental stress regulates invasion front of interactive species. For the 1D case, we tried to see whether any Langevin equation, as predicted by stochastic reaction-diffusion studies, does describe invasion front movements and then obtained corresponding dependencies on environmental stress for cooperative/competitive populations. For the 2D environment, after finding the Langevin equation of invasion front motion, we considered it a growing interface and studied how scaling exponents depend on environmental stress and inter-specific interactions.

Model
Here, we develop an individual-based model in which species live on lattice units. The single-species model in a 1D environment follows these rules: As the initial condition, one cell is located at the first unit. For time evolution, a unit will be selected randomly. Throughout this work, one time step is counted when the number of random selections reaches the number of units defined in the model. If the selected unit does contain a species and there is an empty nearest neighbor (NN) (if the same species occupy both NNs, unit selection will be repeated), then (i) it decides to duplicate into an empty NN and would do so if that unit is occupiable. If the selected NN is not occupiable, the trial number for that NN, n, increases by one (barrier intensity decreases by one). (ii) Independent of the duplication process, the species decides to migrate to an empty NN and do so if the selected NN is occupiable. If the selected NN is not occupiable, the trial number for that unit, n, increases by one (barrier intensity decreases by one). (iii) Any unit would be occupiable after n ≥ N times being selected for migration or duplication where N is environmental barrier intensity (see Fig. 1a). Based on these rules, if a species lives in a unit with two empty and occupiable NNs, each one of them can be occupied by newly created species or through migration. First, the species decides to duplicate and choose one of the NNs to duplicate into it. Since both NNs are empty and occupiable, one of them randomly will be occupied by the species. Then, the initial species decides to migrate. One of the NNs will be selected and if the empty one gets elected, migration happens. If the filled unit is set for migration, the trial fails.
In the two-species model, for simplicity, we consider species to be able to occupy a unit simultaneously. Due to this assumption, they do not compete for space and interaction between species is limited to their mutual effort to degrade the environmental barrier at invasion front 44 . As the initial condition, the first unit is occupied by two species. For time evolution, a unit is selected randomly. If the unit contains one species, it will evolve based on the rules mentioned above (i-iii). If the unit contains both species, one will be selected randomly for duplication and then the other one will be selected for duplication. For migration, again, one of the species will be selected to migrate first. In their interactions with the environment, we count their trials separately as n 1 and n 2 . For an A unit randomly will be selected and duplication and migration trial happen independently. If the selected unit contains a species and two nearest neighbors are empty, the species duplicates into one and can migrate to the other one. The blue arrow shows the direction of upcoming migration to an empty nearest neighbor and the cyan arrow shows upcoming duplication into an empty nearest neighbor. The red arrow shows a failed trial to occupy an empty nearest neighbor. While this attempt has failed, the strength of the barrier has decreased by 1. In a simple singlespecies model, after n = N trials, the unit becomes occupiable. (b) Schematic illustration of the two-species model in one dimension. In a randomly selected unit, each species tries to occupy an empty nearest neighbor. We should have ζ n 1 + n 2 > N or n 1 + ζ n 2 > N for a unit to become occupiable. Migration and duplication happen similar to the single-species model. If the selected unit contains both species, one will be selected randomly for duplication (migration) first and then the other will be selected. However, since we do not include spatial exclusion, selection does not have a relevant role in the majority of cases. www.nature.com/scientificreports/ entirely cooperative scenario, if the number of trials on a unit together exceeds the barrier intensity, n 1 + n 2 > N , the unit becomes occupiable. In a more complex scenario, a unit would be occupiable if we have ζ n 1 + n 2 > N or n 1 + ζ n 2 > N in which ζ is the interaction parameter. When two populations are cooperative (competitive) we have ζ > 0 ( ζ < 0 ). We anticipate a Langevin equation for invasion front, X (see Fig. 1b), and we would try to find out how the diffusion constant and velocity of this interface is related to environmental barrier intensity and interspecific interactions. The 2D version of the model, which is an extension of 1D, will be explained later.

Results
1D case. First, we study the one-dimensional case. We locate a cell at the first unit of a half limited array and let the system evolve based on the above-mentioned rules. We call the occupied unit with the largest distance from the origin as the invasion front (border) location and call its index as X. To find the behavior of invasion front and quantify it, we analyze X, X and X −X where X is the ensemble average of X over different realizations.
Analysis of X −X versus time (Fig. 2a, b) shows that while N affects the magnitude of fluctuations for X −X , the mean squared displacement behaves like a simple random walk and we have: �(X −X) 2 � ∼ t . As a result, we can define a diffusion constant for these fluctuations as �(X −X) 2 � = D f t . Then we studied the dependency of D f on N. It appeared that for large values of N, we have D f ∝ N −γ D with γ D = 2.00 ± 0.05 (Fig. 2c). The averaged velocity of the front position gives us the invasion velocity, V f . Invasion velocity also depends on N as V f ∝ N −γ V with γ V = 1.00 ± 0.05 (Fig. 2c). Such an effect on invasion velocity is expected from an analytical perspective. In a simple 1D invasion model, invasion front velocity, assuming no migration, should be proportional to duplication rate (which is different from that of a fisher's equation, V f ∝ √ RD ). Since duplication happens after N trials, considering the environmental barriers slows down the duplication rate by the factor of N. Respectively, invasion  www.nature.com/scientificreports/ velocity should be slower by the same factor ( N −1 ). Regarding the dependency of the diffusion constant on N, it is enough to look at the definition of D f . As mentioned, invasion velocity decreases by the factor of N.
We now add a second population which does not interfere with the first population except for degrading the barrier in the invasion front. As such, the two populations see each other only on the invasion front. We start the model with two species, located at the first unit and use the same definition for the border, but it does not matter which population has occupied that unit. As mentioned, a unit would be occupiable only if ζ n 1 + n 2 > N or n 1 + ζ n 2 > N . The positive values of ζ show the cooperation between entities and negative values represent the competitive populations. We first try to see how the normalized diffusion constant, N 2 D f , depends on ζ . As Fig. 3a shows, the interaction changes the diffusion constant. Both competitive ( ζ < 0 ) and cooperative ( ζ < 0 ) populations have higher diffusion constant in respect to non-interactive populations ( ζ = 0 ). To see how interaction affects the system response to N, we study the behavior of D f versus N for different interactions (Fig. 3b). Interestingly, the magnitude of diffusion constant depends on interactions, but its behavior versus N, exhibited in value of γ D , depends on interactions. As such, interaction leads to higher diffusion constant with smaller γ D .
We also studied the effect of interactions on invasion velocity, V f . As Fig. 3c shows, cooperation (competition) increases (decrease) the invasion velocity but the effect also is intensified by N. Analysis of behavior of V f versus N shows that γ V also depends on ζ (Fig. 3d). Finally, as Fig. 3e shows, γ V and γ D differently depend on interaction term, ζ . 1D environment may not seem realistic, yet it has played a central role in advancing our understanding of invasion 6,29 . 1D version of our model reveals that competition changes invasion velocity but, for a wide range of interactions ( −0.6 < ζ < 1 ), invasion properties remain primarily unchanged, suggesting that even competition may not affect the invasion velocity significantly. However, to better understand invasion, we need to study the problem in higher dimensions.
2D case. We studied the two-dimension version of our model as well. The 2D version is essentially a semiinfinite array of sites in one direction (just like the 1D case) and a periodic array of sites in the perpendicular direction, with L sites. The same rules will be applied to the 2D case (see Fig. 4a). As the initial condition, all units in the first row will be occupied by species. Migration and duplication can happen into four NNs around each randomly selected unit in single-species and two-species cases. We analyze two different aspects of invasion in 2D environments: invasion velocity and the geometry of the invasion front. The average location of interface, H, is considered as the location of invasion front by setting H =X in which X stands for the average value of X along the border. Similar to 1D, we anticipate a Langevin equation as dH/dt = V H t + Then we studied how H evolve during time to find the interface velocity, V H . As Fig. 4d shows, V H depends on N as V f did.
Due to importance of the geometry of invasion front, we study of the morphology of invasion front through dynamic scaling analysis. For such analysis, we need to calculate surface's width as W 2 = 1 L L i (X i − H) 2 where X i is the invasion front at point i. For variety of surfaces that follow scaling, one has W ≈ L α f (t/L z ) where f(u) is a scaling function such that, f (u) ∝ u β if u ≪ 1 , and f (u) ≈ constant for u ≫ 1 , so that for a fixed L, W ∝ t β . α and β are, respectively, the surface roughness and growth exponents, and z = α/β is the dynamic exponent 45,46 .
We first obtain β for different populations and as Fig. 5a shows, it decreases with N for non-competitive populations. However, for competitive populations, β increases by N. To study the morphology of interface at steady state, we choose the saturated value of interface width, W sat , which is the value W approaches in long times., and analyze its behavior versus L, N and ζ . We tried to determine whether the interface follows any scaling behavior and then see how ζ affects the corresponding exponents. As Fig. 5b shows, W sat for ζ = −1 is larger than other cases, which indicates that invasion front of competitive populations is rougher. Later we studied the behavior of W sat versus L and found the corresponding exponent, α L (or simply α which is the roughness exponent). As Fig. 5c shows, for noncompetitive populations ( ζ ≥ 0 ) α L = 0.70 ± 0.02 but for competitive ones ( ζ = −1 ) we have α L = 0.98 . Finally, we studied the effect of N on W sat . Interestingly, as Fig. 5d shows, competition ( ζ = −1 ) not only leads to higher interface roughness, the associated roughness also is less sensitive to N. While for ζ = −1 we have α N = 0.41 , for ζ ≥ 0 we have α N = 0.61 ± 0.02.
These results reveal that invasion front morphology and its velocity and fluctuations depend on both interspecific interactions and environmental barriers intensity. These results reveal that the geometry of invasion front in competitive populations has significantly different behavior in response to stress. Thus, theoretically, one might be able to estimate clonal interactions by looking at the geometry of the invasion front. In competitive populations, each population cancels the effort of the counterpart population to occupy a new unit. When the intensity of barriers is low, it is more likely for one of the populations to degrade the barrier and occupy a new unit. However, as the intensity increases, such random events become less and less likely, leading to a significant decrease in invasion velocity. It should be noted that adding spatial heterogeneity to barrier resistance can increase the roughness of invasion front 30

Discussion
Understanding tumor invasion through mathematical modeling and in vivo or in vitro studies has significantly increased our understanding of underlying mechanisms. In a now-classic example, it was suggested that duplication rate and diffusion rate of tumor cells determine tumor invasion velocity 47 . Since then, more parameters have been identified and taken to account to understand tumor invasion 48 . Yet, a lot has been remained to be understood about how the interplay between clonal interactions and environmental stress regulate tumor invasion. The fact that geometry of tumor plays a role in patient outcome, puts additional stress on the understanding of the behavior of invasion front. The physical structure of the tumor environment provides a physical barrier against migration and cancer cells need to degrade the ECM to invade 9,11 . As such, we translate the stress to the physical barrier with the intensity of N. Here we considered the physical barrier as a limiting factor that prohibits further growth and cells need to degrade it. We found how interaction significantly regulates invasion velocity. Our results suggest that cooperation plays a crucial role in cells' ability to overcome such a barrier. This conclusion is conceptually in line with other results on the relation between clonal interactions and environmental stress, such as nutrient shortage. We recently showed that once tumor cells individually acquire the ability to induce angiogenesis (angiogenic switch), they may not be able to grow larger until they cooperatively induce further angiogenesis 49 .
The geometry of the invasion front has been used to understand and predict tumor outcome 33 . As a growing interface, scaling analysis has been used to characterize the geometry of invasion front in different studies 30,34,36 . Most of these analyses have concentrated on how environmental features and cellular phenotype and activities such as duplication or dispersal affect the geometry of invasion front 31,34,36 , omitting the direct role of clonal interactions. In this model, two species can simultaneously occupy the same unit. Thus, they do not compete over limited space and their interactions are limited to invasion front. This assumption allows us to highlight the role of interactions on invasion front behavior (adding spatial exclusion makes the problem significantly challenging as we have studied in upcoming work). Here we showed that clonal interactions can single-handedly regulate invasion velocity and the geometry of the invasion front. www.nature.com/scientificreports/ Irregularity of invasion front is associated with of tumor invasive behavior 32,33 , however, the reason behind this association has not been understood yet. Our results here show that local competition between individual cells can lead to irregular invasion front. On the other hand, the relation between poor clinical outcome and clonal diversity is well-established [18][19][20][21][22][23][24] . Thus, irregular geometry is not the cause of adverse clinical outcomes. Instead, both the irregular geometry and adverse outcome are results of clonal diversity and competition.

Summary
Motivated by clonal interactions and environmental barriers that tumor cells experience, we developed a model to study how interspecific interactions and environmental stresses together regulate invasion. In 1D, we found the Langevin equation for invasion front and quantified the dependency of velocity and diffusion constant on the intensity of environmental barriers and the nature of interactions. It turned out that for single-species case, the invasion velocity depends on N as V f ∝ N −γ V with γ V = 1.0 and for the diffusion constant for invasion front, we have D f ∝ N γ D with γ D = 2.0 . Also, competitive populations are more vulnerable to environmental stress and their invasion velocity falls faster in response to N with γ V = 1.80 ± 0.04 . Diffusion constant for interactive populations ( ζ = 0 ) was generally larger and less sensitive to N compared to non-interactive populations ( ζ = 0 or single population model). For the 2D case, the averaged invasion front (H) follows a Langevin equation which depends on N similar to 1D. The geometry of the invasion front exhibits scaling behavior. For ζ = −1 , we found that N increases β . The behavior of W sat versus N, L and ζ was obtained and it turned out that competition not only leads to more rough interfaces, but it also makes those interfaces resistant to environmental stresses. These findings deepen our understanding of the invasion of interactive species and may have applications to understanding tumor clonal interactions during the invasion. www.nature.com/scientificreports/ It should be mentioned that real-world invasions, even those in highly controlled environments, are inherently complex processes. When it comes to cancer, invasion is highly regulated at different levels. Phenotypic plasticity, paracrine interactions with immune cells and fluctuating stresses are only a few relevant examples of huge number of processes that participate in the complex process of invasion. The goal of this work, has not been to deny other biologically relevant factors. Instead, we have tried to concentrate on the role of population level interactions and how they can change invasion front properties by simplifying different aspects.