Papers in Physics, vol. 13, art. 130002 (2021) Received: 17 October 2020, Accepted: 23 February 2021 Edited by: R. Cerbino Reviewed by: F. Giavazzi, Università degli Studi di Milano, Italy Licence: Creative Commons Attribution 4.0 DOI: https://doi.org/10.4279/PIP.130002 www.papersinphysics.org ISSN 1852-4249 Radial percolation reveals that Cancer Stem Cells are trapped in the core of colonies Lucas Barberis1* Using geometrical arguments, it is shown that Cancer Stem Cells (CSCs) must be con- fined inside solid tumors under natural conditions. Aided by an agent-based model and percolation theory, the probability of a CSC being positioned at the border of a colony is estimated. This probability is estimated as a function of the CSC self-renewal probability ps; i.e., the chance that a CSC remains undifferentiated after mitosis. In the most common situations ps is low, and most CSCs produce differentiated cells at a very low rate. The results presented here show that CSCs form a small core in the center of a cancer cell colony; they become quiescent due to the lack of space to proliferate, which stabilizes their population size. This result provides a simple explanation for the CSC niche size, dispensing with the need for quorum sensing or other proposed signaling mechanisms. It also supports the hypothesis that metastases are likely to start at the very beginning of tumor development. I Introduction Cancer Stem Cells (CSCs) are responsible for driv- ing tumor growth due to their ability to make copies of themselves (self-renewal) and differentiate into cells with more specific functions [1]. Differentiated Cancer Cells (DCCs) maintain a limited ability to proliferate, but can only generate cells of their spe- cific lineage [2]. Like the tissues from which they arise, solid tumors are composed of a heterogeneous popula- tion of cells; many properties of normal stem cells are shared by at least a subset of cancer cells [3, 4]. In many tissues, normal stem cells must be able to migrate to different regions of an organ, where they give rise to the specifically differenti- ated cells required by the organism. These features *lbarberis@unc.edu.ar 1 Instituto de F́ısica Enrique Gaviola (IFEG) and Facul- tad de Matemática, Astronomı́a, F́ısica y Computación, CONICET and UNC, X5000HUE Córdoba, Argentina. are reminiscent of invasion and immortality, two hallmark properties of cancer cells [5, 6]. Resis- tance to chemo/radio therapies gives the CSCs a high chance of survival and of forming new tumors, even after treatment [7]. Hence, based on the con- cept that destroying or incapacitating CSCs would be an efficient method of cancer containment and control, new therapeutic paradigms are the focus of current research. A tumorsphere assay is an experimental biolog- ical model used for the study of CSC features. A tumorsphere is a clonal aggregate of cancer cells grown in vitro from a single cell. It has been shown experimentally that DCCs can not generate a tu- morsphere because they are unable to form com- pact long-term aggregates. As a consequence, the current experimental convention for the definition of CSCs, from a functional point of view, is based on the capacity of a single cell, the seed, to grow a more or less spherical aggregate in a gel suspen- sion. This is why CSCs are said to “drive” tumor progression. 130002-1 Papers in Physics, vol. 13, art. 130002 (2021) / L. Barberis In an experimental assay, where the time evolu- tion of the number of cells in a tumorsphere is mea- sured, it is possible to determine a proliferation rate r consistent with the Population Doubling Time (PDT) of the total population. Unfortunately, this quantity cannot discriminate between the growth rates of CSCs and their differentiated counterparts. Furthermore, measuring the PDT of the CSCs is experimentally very complex [8]. Understanding r is crucial for mathematical modeling in a sys- tem where the offspring might belong to a different population than their parents. Indeed, CSC du- plication could have three possible outcomes: stem cell replication (self-renewal), asymmetric differen- tiation and symmetric differentiation. To model such a feature mathematically, we can assume that the three outcomes will occur with probabilities ps, pa and pd, respectively. From the point of view of the populations, the last two possibilities corre- spond to the birth of a new DCC, the parent cell either keeping (in the asymmetric duplication) or losing (in the symmetric case) its stemness. For example, writing the corresponding population dy- namics differential equations, Beńıtez et al. showed that a CSC will give birth to another CSC at a rate rps which is estimated by fitting their mathemati- cal model to the experimental data [9,10]. Since the probability of self-replication ps seems to be small in homeostasis and in most common culture media, then rps will also be small, leading to quiescence, a main feature of CSCs. Nevertheless, CSCs can be experimentally forced to abandon their quiescent state by using specific growth factors that inhibit differentiation [11, 12] or by restricting oxygen con- centration, as discussed below. In these situations, ps is close to one and tumorspheres will contain a high fraction of CSCs. Also, pa is usually small, although its value can increase under abnormal sit- uations such as injury or disease. Metastasis, the invasion process by which cancer cells leave a tumor to form a new colony at another location, is an intriguing feature of cancer disease. Because CSCs are the seeds of tumors, it is ac- cepted that metastatic tumors must come from a pre-existent primary tumor and might be located near its surface in order to detach, migrate and fi- nally invade another site in the organism. The distribution of CSCs in a primary tumor is key to the understanding of metastasis, cell prolifer- ation and drug resistance. Curiously, according to the three-concentric-layer model proposed by Per- sano et al., CSCs seem to be located in the inner core of glioblastomas. On the other hand, cells on the periphery of the tumor show a more differen- tiated phenotype that is highly sensitive to Temo- zolomide, a drug for cancer treatment [13]. Inter- estingly, these authors demonstrate that CSCs are more proliferative under hypoxic conditions. Fur- thermore, Li et al. report that hypoxia plays an important role in the de-differentiation of cells [14]. These results indicate that a hypoxic environment will increase the numbers of CSCs. The aim of this work is to simulate colony growth by means of an Agent-Based Model (ABM) that mimics basic features of CSC proliferation, with emphasis on its geometrical properties. Multi- ple realizations allow us to estimate the fraction of CSCs situated on the periphery of the colony, showing that it is quite small for large, long-lived colonies. We also use some elements of percolation theory to help interpret and quantify the simulation results. We then report a transition to percolation which depends on the self-renewal probability, ps, of the CSC population. Finally, we conclude that that our results lead to a simple explanation for CSC niche size, and support the hypothesis that the metastatic process must start at the very be- ginning of tumor development. II Simulation of two-dimensional colonies Experimental monoclonal colonies start with a CSC in a suitable culture medium. This cell and its daughters duplicate at a more or less constant rate, forming a colony with an approximately circular shape. The reader interested in examples may re- fer to [15–17] for further information. In this work we use the same principle, describing the cells us- ing circle-shaped mathematical objects and refining rules for their duplication and movement. In later studies we will requirecells not to be rigid spheres. Each cell in the current simulations has a central circular impenetrable region, the core, and around the core there is an area, the corona, where overlap- ping is allowed. In this way, cells are able to overlap as far as the border of their coronas and reach the border of the core of other, neighboring cells. The figures presented in this work has a 95 % overlap 130002-2 Papers in Physics, vol. 13, art. 130002 (2021) / L. Barberis (a) (b) Figure 1: Two realizations of the simulation for ps = 0.5 after 15 days. (a) There are no CSCs at the periphery; all of them, colored in pink, became quiescent. (b) The active CSCs at the border are colored in red. Note that some branches failed to percolate. DCCs are depicted in cyan (quiescent) and blue (active). The seed is rep- resented in yellow. on the cell radius. However, we will show that the results presented here do not depend on this param- eter. Also, each cell belongs to a class that defines its behavior: active CSCs, active DCCs, quiescent CSCs and quiescent DCCs. As a simplification of the model, we assume that the probability of asym- metric duplication of CSCs is zero, pa = 0. In this way ps is the only control parameter that allows direct capture of the model’s main features. We start the simulation by seeding an active CSC and asking it to duplicate. Furthermore, at each time step we ask all active cells to attempt to du- plicate, according to the following rules. First a new cell is created in a random position at the side of an active cell. If there are no other cells in this location it remains there, otherwise a new location beside the active cell is randomly determined. After 500 attempts, if the new cell has still not found an empty spot to occupy, it is destroyed and the active cell changes to its respective quiescent class. Con- versely, after successful duplication, if the active cell is a CSC it self-renews with a probability ps, which means that the new cell becomes an active CSC, otherwise it changes to the active DCC class. If the new, active cells belong to different classes, there is a probability of 1/2 that the active CSC will exchange positions with the new DCC. Natu- rally, an active DCC can only create another active DCC. Once all the active cells have been asked to du- plicate, independently of the success of the attempt (a) ps = 0.2 (b) ps = 0.95 Figure 2: Two realizations of the simulation at low (a) and high (b) self-replication rates after 15 days. In (a) the CSCs are quickly surrounded by DCCs, becoming quiescent (pink). In (b) the high number of CSCs allows percolation to most of the colony perimeter. The seed is colored in yellow. we say that a one-day-long time step has been per- formed. This implies, without loss of generality, that r ' 1 day−1, a reasonable growth rate that will aid our intuition. Further details on the struc- ture of the simulation process are provided in the Appendix, including a flow chart in Fig. 7. We have also provided videos in the Supplementary In- formation. In each video a representative value of ps was set and ten realizations were recorded, illus- trating how they were carried out. Typical outcomes of the described process are shown in Fig. 1 and video v1. In these examples we set ps = 0.5 and started with a CSC seed depicted in yellow. After running the simulation for a period equivalent to two weeks (15 time steps), a relatively long period for a biological experiment, we obtained two possible outcomes: (a) a few quiescent CSCs were trapped in the center of the colony or (b) there were active CSCs at the border of the colony. In- deed, we defined the border of the colony as the rim formed by the set of active cells. To better track the subpopulations present in the colony, the ac- tive cells are colored in red for CSC and in blue for DCC. Also, the quiescent cells are colored in pink for CSC and cyan for DCC. Remember that the yellow dot is not always at the center of the colony because the seed may exchange its pposition with a new DCC. In panel (b) we recognize a path, paved by CSCs, that joins the center of the colony with its border. Also, some frustrated branches appear in this path of CSCs, which die in the quiescent core, indicating great variability. This path resembles 130002-3 https://www.papersinphysics.org/papersinphysics/article/view/641/349 Papers in Physics, vol. 13, art. 130002 (2021) / L. Barberis 0 10 20 30 40 50 t 10 0 10 1 10 2 A c ti v e C S C (a) 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 10 20 30 40 50 t 10 0 10 1 10 2 10 3 10 4 T o ta l C S C (b) 0 10 20 30 40 50 t 0 0.2 0.4 0.6 0.8 1 P (c) 0 10 20 30 40 50 t 10 0 10 1 10 2 10 3 10 4 N (d) Figure 3: Time evolutions of quantities averaged over a thousand realizations for values of ps ∈ [0.2..1.0]. (a) Number of CSCs at the border. (b) The total number of CSCs reaches a constant value that defines the niche. (c) The probability of a CSC being at the border. (d) System size as a function of time. cluster percolation in porous media, lattices or net- works, which led to us to attempt to describe the probability of finding a CSC at the border of the colony by means of percolation theory. To sharpen our intuition as to what percolation means in this system, in Fig. 2 and videos v2 and v3 we present two more examples with different self-replication probabilities, ps. In panel (a) and video v2, ps = 0.2 is so small that the CSCs are quickly surrounded by DCCs, and become quies- cent. This is the most common situation in a regu- lar culture medium where the CSC count is low and constant —over time. On the other hand, in panel (b) and video v3, ps = 0.95 leads to a large CSC population that invades almost the entire system. The addition of stem cell maintenance factors such as EGF or bFGF to the culture medium is an ex- ample of this case. These limiting situations were previously studied both experimentally [11, 18] and mathematically [9, 10], focusing in the first case on technical aspects of the assay, and in the second case on recovery of the CSC fraction. Curiously, in the examples shown in Figs. 1 and 2, the rim formed by active cells is similar to the one reported in the classical multicelular spheroid work of Freyer and Shuterland [19]. In Freyer’s work, the inner portion of the spheroids was formed by dead cells, attributed to a low concentration of oxy- gen/nutrients supplied by means of diffusion. This rim was also formed by active cells and was stud- ied mathematically by several authors using diffu- sion equations [20–22]. It is interesting that in the present monoclonal case the geometry seems to be sufficient to develop a rim, just two cells thick, independently of any diffusion process. 130002-4 https://www.papersinphysics.org/papersinphysics/article/view/641/350 https://www.papersinphysics.org/papersinphysics/article/view/641/351 https://www.papersinphysics.org/papersinphysics/article/view/641/350 https://www.papersinphysics.org/papersinphysics/article/view/641/351 Papers in Physics, vol. 13, art. 130002 (2021) / L. Barberis To study the percolation properties of the system statistically, we carried out extensive simulations of colony growth. Each data point reported is the result of averaging over 1000 runs. The time evolu- tion of colony growth was performed for 37 values of ps ∈ [0.2, 1.0]. The results of these measure- ments are summarized in Fig. 3. Panel (a) shows the number of CSCs at the border, which increases to a maximum given by the time when they are overwhelmed by the DCC population. Beyond this point, more and more CSCs become quiescent un- til no more active CSCs remain. This phenomenon disappears for ps barely lower than 1, when almost all cells on the border are CSCs. In panel (b) the time evolution of the total CSC population is de- picted, showing that after a transient the CSC pop- ulation stops growing and remains constant. This fixed number of CSCs is usually called the CSC niche and was mathematically studied in [9]. As expected, the probability of a CSC being at the border rises as ps is increased, as shown in panel (c). The relationship between simulation time and system size is shown in panel (d); due to its geomet- rical nature it is independent of ps, the self-renewal probability. III Percolation theory Percolation theory was used to estimate the prob- ability of finding a CSC at the border of the colony, looking for a purely geometrical feature of its growth. In particular, we assumed that the cells in the colony could be mapped onto the nodes of a triangular network with the seed at its center. Each cell is then connected with its nearest neigh- bors, whose number, in two dimensions, is known not to exceed six. As shown in Fig. 1(b), the CSCs (quiescent and active) form paths that expand ra- dially from the center to the edge of the colony. By inspection, we note that these paths are formed by the connected nodes that have had a CSC at any moment during the culture. The path could have some “holes” that arise because of the possi- bility of exchanging a parent CSC with its daughter DCC. The probability ps thus regulates the num- ber of CSCs that occupy the nodes of the network. Hence, there must be a threshold value pc of ps, be- low which it is impossible to follow a path of nodes starting at the seed and ending in a CSC at the border of the tumor. This value, pc, is called the critical value or percolation threshold and defines a percolation phase transition[23]. Thus, when the border is reached by CSCs, we say that such a path percolates. Formally, a percolating path is mathematically defined as a set of connected nodes that expands in- finitely over an infinite network for a large enough value of ps [24]. In our model, each new node of the percolation path is added either by the self- renewal of a CSC, with probability ps, or by an exchange between a CSC and a DCC, with proba- bility 1 2 ×(1−ps), after a CSC gives birth to a DCC. It is also important to note that several neighboring nodes could already be occupied. Thus, it is useful to define the average empty neighbor number z as follows: Imagine that we perform a random walk starting at the seed’s node along an infinite path where it is forbidden to step back. At each step there are z branches in the network, but there is also a probability ps + 1−ps 2 that a node belongs to the path, so on average only 1 2 (ps + 1)z will be ac- cessible. To continue walking along the path, there must be at least one empty node to walk along, leading to the condition 1 2 (ps + 1)z ≥ 1. Therefore, the critical probability for transition to percolation is pc = 2 z − 1, (1) which depends on the available neighbor number. 1 It is noteworthy that in Eq. (1), for z = 1 we ob- tain pc = 1, the 1D critical occupation for a chain of nodes. For 1 < z < 2 the critical occupation quickly falls to zero. This means that for z ≥ 2 the system will always percolate. On the other hand, if the CSC never exchanges its position with a DCC, we get pc = 1 z , which is again the 1D critical oc- cupation and will asymptotically fall to zero as z increases. Otherwise, if the CSC always moves, ps diverges as expected for a cell that always remains on the outer rim of the growing colony. There is no evidence that a CSC will stay either in the core or on the border of the colony. Simulations car- ried out on both limits of the exchange probability agree with these theoretical limits, but do not pro- vide further information. For this reason, we im- plemented the one-half exchange probability in our 1Remember it is not universal because it depends on the details of the network. 130002-5 Papers in Physics, vol. 13, art. 130002 (2021) / L. Barberis (a) (b) 0 10 20 30 40 50 n 1 1.2 1.4 1.6 1.8 2 z Figure 4: (a) Cell overlapping does not change the max- imum neighbor number of six; thus, the underlying lat- tice is topologically triangular (yellow lines). (b) The average neighbor number versus layer number quickly drops to one. At t = 50 days, n ' 50, N ' 4000 and z = 1.02, which means ps = 0.92. simulations, using the “a priori equal probability” principle. By construction, there are no closed loops in the network, thus deduction of the critical probability was similar to that for a Bethe lattice or a tree graph. However, unlike regular lattices/graphs, in our model the neighbor number of a node is not the same for all nodes due the randomness intro- duced by space searching and the probability of exchange. For this reason, we defined an average neighbor number. In Fig. 4(a) each cell is depicted by two concentric circles; the dark one represents the rigid core and the lighter represents the corona where overlapping is allowed. Note that the maxi- mum neighbor number is six, regardless of the ex- tent of the corona. In simulations a new cell must be in contact with other cells, and overlapping will modify the density of the colony. As mentioned, this parameter will be useful for further study of diffusion effects that are irrelevant to the present work. Also, the random process of searching space for proliferation will slightly modify the geometry of the underlying network, but not its topology, which will be the same as that of the triangular lattice except for a few defects caused by a neigh- bor number lower than six. Thus density, which is relevant when studying the diffusion of nutrients, oxygen or proteins for signaling, does not play any role in the current percolation problem. To roughly estimate z we built up a tree graph Figure 5: Growth of the colony on a triangular lattice, starting from the center (black dot). Each layer is de- picted in a different color and contains three nodes with z = 3 and the remainder with z = 1. over a triangular lattice substrate. A possible con- nection pattern of nodes in such a network is de- picted in Fig. 5, which was built layer by layer with each layer represented in a different color. Layer n = 0 is the central black dot that represents the seed. In this particular example, in each layer we depict three nodes with z = 3 (darker dots), leav- ing the remaining nodes with z = 1 (lighter dots). It is clear that for any layer there are 6 ×n nodes for n > 0; then to connect two layers without loops we need on average z = n+1 n . Because 1 < z < 2 quickly decays with n, c.f. Fig. 4(b), we do expect that the critical probability for percolation, given by Eq. (1), will shift to pc = n n+1 → 1 as n → ∞. The size of the network will be N = 6 ∑n i=1 i + 1 leading to the limit pc −−−−→ N→∞ 1. Thus, as the colony increases in size, the percolation threshold shifts towards 1. The main consequence is that there is no chance to perform finite size scaling to detect the critical transition point through simula- tions, as in the case of 1D chains. As mentioned before, in Figs. 1 and 2, the pro- liferative cells are distributed in the two outermost layers. This feature of the colony occurs because the cells in the simulation are requested to prolif- erate in random order. This is also the cause of 130002-6 Papers in Physics, vol. 13, art. 130002 (2021) / L. Barberis 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 p s 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 P t [days] - N [cells] 20 - 500 40 - 2389 60 - 5797 80 - 10669 100 - 17090 Bethe 0 1 2 10 4 0.2 0.4 0.6 0.8 1 p c N simulation theory Figure 6: Probability of percolation as a function of ps for different colony sizes. Dots correspond to sim- ulations and lines are theoretical fittings to estimate the percolation threshold (see text). The gray line is the theoretical Bethe result. Inset: Percolation threshold pc at different sizes measured by simulations (blue). Our simple theoretical approach overestimates this threshold (red). the observed circular shape of the colony. Because the vertices of the hexagonal array shown in Fig. 5 have a lower probability of being chosen for prolifer- ation, as system size increase s cells on their edges will proliferate first, as at the beginning of each time step they have more room to do so. Thus, the estimated number of cells in each layer is under- estimated, as is their neighbor number. As a con- sequence, we expect that the actual value of pc as a function of size or time will be much lower than that predicted by Eq. (1) and our estimation of z. A precise estimation of z will require statistical treatment that is beyond the scope of this article. The result implies that proper (mathematical) percolation transition will occur at ps = 1, as in 1D percolation. This also implies that the only way for a CSC to be at the border, for an arbitrarily large colony, is when it is forced not to differentiate. Nevertheless, monoclonal colonies cannot be main- tained and grown forever. A 50-day experiment is very long, and very uncommon in the literature. We thus focus on the probability of a CSC being at the border, P∞, whose behavior, as a function of ps, is depicted in Fig. 6 for several times/sizes. The relationship between time and size is bijective, thus we report both values in the legend. In this graph the dots are the result of the simulations pre- sented in the previous section, for which the fre- quency of CSCs at the border of the colony was averaged over a thousand runs. To obtain suitable fitting of these data as continuous functions of ps, we used erf functions following Yonezawa’s classi- cal work [23]. The inflection point of the erf curves coincides with the percolation threshold, giving a good estimate of pc. In addition, the theoreti- cal result for the percolation parameter in a Bethe lattice of order 3 is depicted in gray for compar- ative purposes. Note that P∞ does not jump as in Bethe percolation; the transition is continuous, as known for several regular lattices, including the triangular one. Gonzalez et al. [25] studied the transitions of several forms of percolation in trian- gular lattices, reporting their critical probabilities using finite size scaling analysis. They found uni- versal features with values between 0.5 and 0.8 for the percolation threshold, depending on the prob- lem. In contrast, but as expected, as our system size increases the fitted curves become steeper and steeper, and their inflection points predict a shift of ps towards 1, as depicted by the blue dots in the inset of Fig. 6. Comparison with our theoret- ical result gives an expected overestimation of the percolation threshold – the red line, as explained before. We therefore hope to find that a CSC will be at the border of the colony in half the runs when we simulate a three-week experiment. IV Conclusions Metastasis is an intriguing feature of cancer inva- sion. Most research in this field is guided by the bi- ological point of view, even though there are many mathematical models that attempt to describe it [26–29]. In the examples of growth given in Fig. 1(a) and Fig. 2(a), it becomes clear that under normal cul- ture conditions CSCs remain in the inner core of spheroids. This is also deduced from the low P∞ probability of finding a CSC at the border of the colony, as shown in Fig. 6. If this could be observed in the laboratory, the experiments that reveal a CSC preference for hypoxic environments could be 130002-7 Papers in Physics, vol. 13, art. 130002 (2021) / L. Barberis easily explained by this fact. It has been shown experimentally that hypoxia maintains the undif- ferentiated state of primary glioma cells, slowing down the growth of glioma cells which are in a rel- atively quiescent stage, increasing colony-forming efficiency and the migration of glioma cells, and ele- vating the expression of of stem cell markers. How- ever, the expression of markers for stem cell dif- ferentiation was reduced after hypoxia treatment [13, 14]. It is also known that the inner core of spheroids and tumors have a lack of oxygen and nutrients, which are first consumed by the outer layers of the tumor [30–33]. In this context, CSCs constitute a phenotype that has evolved to survive under hypoxic conditions in order to drive tumor growth. On the other hand, metastasis requires CSCs to overcome three major barriers: the probability of an active CSC being present at the border of the colony, the probability of it detaching from the tu- mor and the probability of it finding a suitable place to proliferate. The results of the present work es- tablish that the first of these probabilities is very small, even for relatively small tumors, supporting the rare occurrence of metastasis. For large self-replication rates, possibly gener- ated by environmental conditions, many active CSCs would become candidates for detachment. At present, we do not know of any report of a high number of CSCs in tumors in vivo or in transplanted xenographs under normal physiologi- cal conditions that experimentally support a large6 self-replication rate as part of the metastasis mech- anism. An intriguing possibility is that metastasis begins in the early stages of tumor progression [34–36]. This is deduced from Fig. 3(a) where the number of active CSCs shows a peak for low ps and short times. If this is the actual situation, a CSC must detach from the primary tumor in its first week of life, when the probability of being at the border is close to 1, even for low values of self-replication. Note that at this time the tumor consists of no more than a dozen cells. This fact leads us to hy- pothezise that a primary tumor is indeed the metas- tasis of some small young tumors. According to this hypothesis the primary tumor is in fact the first colony able to fix and grow, while metastatic tumors come from later CSCs that take more time to attach successfully to a new location. The de- tached cells must survive while competing for nutri- ents and space with the primary spheroids, which at this stage are ruthless adversaries [37]. Our two-dimensional model could be extended to three dimensions; this would be more realistic but also computationally more expensive. Preliminary results are qualitatively similar to those reported here, butso far we lack the statistics to perform a quantitative comparison. Another feature, ob- tained by changing geometrical constraints in our code, are simulations of non-solid tumors such as haematologic neoplasm; preliminary results show that active CSCs are present in a larger proportion there than in solid tumor spheroids. As a conse- quence, there is a higher probability of CSCs de- taching and developing metastasis in neoplasms. In summary, percolation provides a way of de- veloping a geometrical theory to support or com- plement signaling pathways, quorum sensing and other tools frequently used to study metasta- sis. Further experimental research will elucidate whether CSCs are the survivors with greatest fit- ness, as suggested by our present results. Acknowledgements - This work was supported by SECyT-UNC (project 05/B457) and CONICET (PIP 11220150100644), Argentina. The author is grateful to Dr. C. A. Condat and Dr. L. Vellón for fruitful discussions. Also, he thanks Dr. Fabio Giavazzi for suggesting a better approach for Eq. (1). Appendix: Simulation algorithm. The simulations were carried out in an object- oriented paradigm. Each cell belongs to one of the following classes: Active CSC, active DCC, quies- cent CSC or quiescent CSC. First we created an active CSC object at the center of a square space large enough to contain the final colony. Then for each time step we requested, in random order, that all objects belonging to an active class perform the duplication procedure. In the left-hand column of Fig. 7 we depict the flow chart of the duplication procedure, whose steps are illustrated in the right- hand column. 130002-8 Papers in Physics, vol. 13, art. 130002 (2021) / L. Barberis Move towards Is the Draw a random Is thereDid this 500 times? A random direction is defined The new cell moves until there is no core superposition A parent DCC always creates a DCC If a parent CSC creates a CSC When a CSC creates a DCC and doesn't change position When a CSC creates a DCC and changes position Set new as CSC Set parent as DCC Active DCCAny cell Active CSC Quiescent DCC Quiescent CSC Destroy the new cell Set the parent as quiescent Set new as CSC Set new as DCC Create a new cell Looking for space subroutine A new cell is created in the position of the parent cell Move back Draw a random Set new as DCC Figure 7: Left: Flow chart for the proliferation process of each cell. Right: Illustration of the stages and outcomes of the duplication process. The duplication process starts when an active cell, the parent, creates a copy of itself, the new cell, in the same position. Then a random angle 0 ≤ rα < 2π is drawn from a uniform distribution and the new cell will advance a certain distance in the direction defined by this angle, such that its border does not overlap with its parent’s core. This distance was determined as 95% of the cell diame- ter. If in its new position the border of the new cell overlaps with the core of another cell, it moves back to its initial position. Another random direction is then defined, repeating the movement sequence. The new cell will move forward and back, randomly changing its direction until there is no overlapping; thus we say that an empty spot is found. After 500 failed attempts to find an empty spot, the new cell is destroyed and the parent cell changes to its re- spective quiescent class. Conversely, if the new cell 130002-9 Papers in Physics, vol. 13, art. 130002 (2021) / L. Barberis finds an empty spot, and if its parent is a DCC, it sets its own class as an active DCC. On the other hand, if the parent is a CSC, a new random num- ber 0 ≤ rd < 1 is drawn to be compared with ps. If rd < ps the new cell will set its class to active CSC, otherwise the new cell class will be set to ac- tive DCC. In the latter case a new random num- ber 1 ≤ rm < 1 is drawn, and if the comparison rm < 1/2 becomes true, the cells will exchange their positions, leaving the CSC in the new spot and the DCC in the position of its parent. Note that quiescent cells will never be requested to do anything; they do not move or duplicate in these simulations. Their only role is to restrict the movement of active cells. As a consequence, the code runs very rapidly, even for large colonies. [1] E Batlle, H Clevers, Cancer stem cells revis- ited, Nat. Med. 23, 1124 (2017). [2] C A M La Porta, S Zapperi, J P Sethna, Senescent cells in growing tumors: Popula- tion dynamics and cancer stem cells, PLoS Comput. Biol. 8, e1002316 (2012). [3] N A Lobo, Y Shimono, D Qian, M F Clarke, The biology of cancer stem cells, Annu. Rev. Cell Dev. Biol. 23, 675 (2007). [4] J Stingl, C Caldas, Molecular heterogeneity of breast carcinomas and the cancer stem cell hypothesis, Nat. Rev. Cancer 7, 791 (2007). [5] Y Shimono, M Zabala, R W Cho, N Lobo, P Dalerba, D Qian, M Diehn, H Liu, S P Panula, E Chiao, F M Dirbas, G Somlo, R A Reijo Pera, K Lao, M F Clarke, Downregula- tion of miRNA-200c links breast cancer stem cells with normal stem cells, Cell 138, 592 (2009). [6] D Hanahan, R A Weinberg, The hallmarks of cancer review, Cell 100, 57 (2000). [7] P Jagust, B de Luxán-Delgado, B Parejo- Alonso, P Sancho, Metabolism-based thera- peutic strategies targeting cancer stem cells, Front. Pharmacol. 10, 203 (2019). [8] A Waisman, F Sevlever, M Eĺıas Costa, M S Cosentino, S G Miriuka, A C Ventura, A S Guberman, Cell cycle dynamics of mouse embryonic stem cells in the ground state and during transition to formative pluripotency, Sci. Rep 9, 8051 (2019). [9] L Beńıtez, L Barberis, C A Condat, Mod- eling tumorspheres reveals cancer stem cell niche building and plasticity, Physica A 533, 121906 (2019). [10] L Beńıtez, L Barberis, C A Condat, Un- derstanding the influence of substrate when growing tumorspheres, BMC Cancer (In press) (2021). [11] J Wang, X Liu, Z Jiang, L Li, Z Cui, Y Gao, D Kong, X Liu, A novel method to limit breast cancer stem cells in states of quiescence, pro- liferation or differentiation: Use of gel stress in combination with stem cell growth factors, Oncol. Lett. 12, 1355 (2016). [12] Y C Chen, P N Ingram, S Fouladdel, S P Mc- Dermott, E Azizi, M S Wicha, E Yoon, High- throughput single-cell derived sphere forma- tion for cancer stem-like cell identification and analysis, Sci. Rep. 6, 27301 (2016). [13] L Persano, E Rampazzo, A Della Puppa, F Pistollato, G Basso, The three-layer concentric model of glioblastoma: Cancer stem cells, microenvironmental regulation, and therapeutic implications, TheScientific- WorldJo. 11, 1829 (2011). [14] P Li, C Zhou, L Xu, H Xiao, Hypoxia enhances stemness of cancer stem cells in glioblastoma: An in vitro study, Int. J. Med. Sci. 10, 399 (2013). [15] C Zhang, Y Tian, F Song, C Fu, B Han, Y Wang, Salinomycin inhibits the growth of colorectal carcinoma by targeting tumor stem cells, Oncol. Rep. 34, 2469 (2015). [16] S Shankar, D Nall, S-N Tang, D Meeker, J Passarini, J SharmaR, K Srivastava, Resveratrol inhibits pancreatic cancer stem cell characteristics in human and KrasG12D transgenic mice by inhibit- ing pluripotency maintaining factors and epithelial-mesenchymal transition, PLoS ONE 6, e16530 (2011). 130002-10 https://doi.org/10.1038/nm.4409 https://doi.org/10.1371/journal.pcbi.1002316 https://doi.org/10.1371/journal.pcbi.1002316 https://doi.org/10.1146/annurev.cellbio.22.010305.104154 https://doi.org/10.1146/annurev.cellbio.22.010305.104154 https://doi.org/10.1038/nrc2212 https://doi.org/10.1016/j.cell.2009.07.011 https://doi.org/10.1016/j.cell.2009.07.011 https://doi.org/10.1016/S0092-8674(00)81683-9 https://doi.org/10.3389/fphar.2019.00203 https://doi.org/10.1038/s41598-019-44537-0 https://doi.org/10.1016/j.physa.2019.121906 https://doi.org/10.1016/j.physa.2019.121906 https://doi.org/10.21203/rs.3.rs-67713/v1 https://doi.org/10.21203/rs.3.rs-67713/v1 https://doi.org/10.3892/ol.2016.4757 https://doi.org/10.1038/srep27301 https://doi.org/10.1100/2011/736480 https://doi.org/10.1100/2011/736480 doi:10.7150/ijms.5407 doi:10.7150/ijms.5407 https://doi.org/10.3892/or.2015.4253 https://doi.org/10.1371/journal.pone.0016530 https://doi.org/10.1371/journal.pone.0016530 Papers in Physics, vol. 13, art. 130002 (2021) / L. Barberis [17] A Schneider, D Spitkovsky, P Riess, M Mol- canyi, N Kamisetti, M Maegele, J Hescheler, U Schaefer, “The good into the pot, the bad into the crop!” — A new technology to free stem cells from feeder cells, PLoS ONE 3, e3788 (2008). [18] A Chen, L Wang, S Liu, Y Wang, Y Liu, M Wang, Attraction and compaction of mi- gratory breast cancer cells by bone matrix pro- teins through tumor-osteocyte interactions, Sci. Rep. 8, 5420 (2018). [19] J P Freyer, R M Sutherland, Regulation of growth saturation and development of necro- sis in EMT6/Ro multicellular spheroids by the glucose and oxygen supply, Cancer Res. 46, 3504 (1986). [20] C A Condat, S A Menchón, Ontoge- netic growth of multicellular tumor spheroids, Physica A 371, 76 (2006). [21] S A Menchón, C A Condat, Cancer growth: Predictions of a realistic model, Phys. Rev. E 78, 022901 (2008). [22] L Barberis, C A Condat, Describing interac- tive growth using vector universalities, Ecol. Model. 227, 56 (2012). [23] F Yonezawa, S Sakamoto, M Hori, Percola- tion in two-dimensional lattices. I. A tech- nique for the estimation of thresholds, Phys. Rev. B 40, 636 (1989). [24] K Christensen, N R Moloney, Complexity and Criticality, Imperial College Press, London, UK (2005). [25] M I González, P Centres, W Lebrecht, A J Ramirez-Pastor, F Nieto, Site-bond percola- tion on triangular lattices: Monte Carlo sim- ulation and analytical approach, Physica A 392, 6330 (2013). [26] K Iwata, K Kawasaki, N Shigesada, A dy- namical model for the growth and size dis- tribution of multiple metastatic tumors, J. Theor. Biol. 203, 177 (2000). [27] S A Menchón, C A Condat, Modeling tumor cell shedding, Eur. Biophys. J. 38, 479 (2009). [28] J B McGillen, E A Gaffney, N K Martin, P K Maini, A general reaction-diffusion model of acidity in cancer invasion, J. of Math. Biol. 68, 1199 (2014). [29] A Rhodes, T Hillen, A mathematical model for the immune-mediated theory of metasta- sis, J. Theor. Biol. 482, 109999 (2019). [30] A I Liapis, G G Lipscomb, O K Crosser, E Tsiroyianni-Liapis, A model of oxygen dif- fusion in absorbing tissue, Math. Modell. 3, 83 (1982). [31] P J Robinson, S I Rapoport, Model for drug uptake by brain tumors: Effects of osmotic treatment and of diffusion in brain, J. Cerebr. Blood F. Met. 10, 153 (1990). [32] Y Jiang, J Pjesivac-Grbovic, C Cantrell, J P Freyer, A multiscale model for avascular tu- mor growth, Biophys. J. 89, 3884 (2005). [33] J A Bull, F Mech, T Quaiser, S L Waters, H M Byrne, Mathematical modelling reveals cellular dynamics within tumour spheroids, PLoS Comput. Biol. 16, e1007961 (2020). [34] J W Gray, Evidence emerges for early metas- tasis and parallel evolution of primary and metastatic tumors, Cancer Cell 4, 4 (2003). [35] H Hosseini, M M Obradović, M Hoffmann, K L Harper, M S Sosa, M Werner-Klein, L K Nanduri, C Werno, C Ehrl, M Ma- neck, N Patwary, G Haunschild, M Gužvić, C Reimelt, M Grauvogl, N Eichner, F We- ber, A D Hartkopf, F A Taran, S Y Brucker, T Fehm, B Rack, S Buchholz, R Spang, G Meister, J A Aguirre-Ghiso, C A Klein, Early dissemination seeds metastasis in breast cancer, Nature 540, 552 (2016). [36] N Linde, M Casanova-Acebes, M S Sosa, A Mortha, A Rahman, E Farias, K Harper, E Tardio, I Reyes Torres, J Jones, J Con- deelis, M Merad, J A Aguirre-Ghiso, Macrophages orchestrate breast cancer early dissemination and metastasis, Nature Com- mun. 9, 1 (2018). [37] L Barberis, M A Pasquale, C A Condat, Joint fitting reveals hidden interactions in tumor growth, J. Theor. Biol. 365, 420 (2015). 130002-11 https://doi.org/10.1371/journal.pone.0003788 https://doi.org/10.1371/journal.pone.0003788 https://doi.org/10.1038/s41598-018-23833-1 https://cancerres.aacrjournals.org/content/46/7/3504.short https://cancerres.aacrjournals.org/content/46/7/3504.short https://doi.org/10.1016/j.physa.2006.04.082 https://doi.org/10.1103/PhysRevE.78.022901 https://doi.org/10.1103/PhysRevE.78.022901 https://doi.org/10.1016/j.ecolmodel.2011.12.011 https://doi.org/10.1016/j.ecolmodel.2011.12.011 https://doi.org/10.1103/PhysRevB.40.636 https://doi.org/10.1103/PhysRevB.40.636 https://doi.org/10.1142/p365 https://doi.org/10.1142/p365 https://doi.org/10.1016/j.physa.2013.09.001 https://doi.org/10.1016/j.physa.2013.09.001 https://doi.org/10.1006/jtbi.2000.1075 https://doi.org/10.1006/jtbi.2000.1075 https://doi.org/10.1007/s00249-008-0398-5 https://doi.org/10.1007/s00285-013-0665-7 https://doi.org/10.1007/s00285-013-0665-7 https://doi.org/10.1016/j.jtbi.2019.109999 https://doi.org/10.1016/0270-0255(82)90014-8 https://doi.org/10.1016/0270-0255(82)90014-8 https://doi.org/10.1038/jcbfm.1990.30 https://doi.org/10.1038/jcbfm.1990.30 https://doi.org/10.1529/biophysj.105.060640 https://doi.org/10.1371/journal.pcbi.1007961 https://doi.org/10.1016/S1535-6108(03)00167-3 https://doi.org/10.1038/nature20785 https://doi.org/10.1038/s41467-017-02481-5 https://doi.org/10.1038/s41467-017-02481-5 https://doi.org/10.1016/j.jtbi.2014.10.038 Introduction Simulation of two-dimensional colonies Percolation theory Conclusions