MicroRNA-222 Regulates Melanoma Plasticity

Melanoma is one of the most aggressive and highly resistant tumors. Cell plasticity in melanoma is one of the main culprits behind its metastatic capabilities. The detailed molecular mechanisms controlling melanoma plasticity are still not completely understood. Here we combine mathematical models of phenotypic switching with experiments on IgR39 human melanoma cells to identify possible key targets to impair phenotypic switching. Our mathematical model shows that a cancer stem cell subpopulation within the tumor prevents phenotypic switching of the other cancer cells. Experiments reveal that hsa-mir-222 is a key factor enabling this process. Our results shed new light on melanoma plasticity, providing a potential target and guidance for therapeutic studies.


Introduction
Plasticity of tumor cells is not only a key ingredient of tumor growth but also a critical step for developing a successful anti-tumor therapy. Recently, our group showed that a network of miRNAs modulates phenotypic switching in human melanoma through Wnt signalling and epithelial-mesenchymal transition (EMT) in human melanoma cells, controlling the total number of cancer stem cells (CSCs) in the tumor [1]. Among the identified miRNAs, hsa-mir-222-5p plays a critical role [1]. The ability of CSCs to keep their number constant resembles a critical biological feature of stem cells (SCs). Since SCs live in a confined space (the niche), they need to ensure that their number is kept constant [2]. If for some reason a drastic depletion of SCs occurs, other more differentiated cells can revert to the SC state, replenishing the population. During this process, the CSCs subpopulation overshoots to a concentration higher than the original one and then returns to the initial level. The classical cell hierarchy within the tissue considers SCs at the top of a lineage and their differentiation as an irreversible process. In recent years, this paradigm has been challenged by numerous studies demonstrating that progenitors and differentiated cells in different tissues, from intestine and gastric cells to skin and lung, display considerable plasticity under various conditions [3][4][5][6][7][8]. Furthermore, it has been shown that niche-derived signals contribute to this plasticity [9][10][11][12]. Recent work indeed showed that microenvironmental signals [13] and the density of the extracellular matrix [14] modulate the spatial patterning of CSCs. One of the important consequences of tumor plasticity lies in its impact on drug treatment [15]. Interestingly, a recent paper showed that after continuous drug treatment, slow-cycling melanoma cells enriched within the tumor show a SC-like phenotype [16].
In this paper, we investigated the regulatory mechanisms controlling phenotypic switching in human melanoma IgR39 cells, by combining mathematical models and experiments. We first developed a mathematical model to illustrate how phenotypic switching of cancer cells (CCs) into CSCs can be controlled by the release of switch-inhibiting factors from CSCs. Previous work from our group showed that CSCs in melanoma can be identified using the CXCR6 marker [1,17]. We have also already demonstrated that a stochastic mechanism displaying an environment-independent switching rate would not account for the observed overshoot and we proposed a model which introduces a switching rate that depends on an activating miRNA, whose dynamics are controlled by the exponential of the fraction of CSCs, in agreement with experiments [1]. Even if this model could reproduce experimental observations, it was difficult to translate it into a specific molecular mechanism. In the literature, there are other models that use rate equations to describe the growth of the tumor mass, including phenotypic switching [18][19][20][21][22]. For example, in [19] a key ingredient to describe the de-differentiation of epithelial cells to the more aggressive mesenchymal phenotype is a time delay of some days, resulting in delayed differential equations.
Here we describe an alternative model to control phenotypic switching that can easily be interpreted in terms of molecular interactions in agreement with previous experimental data. Unlike other models, we suggest that phenotypic switching is a spontaneous tendency of all cancer cells, and that, under standard conditions of tumor growth, switching is inhibited by a specific molecular species produced by CSCs. We then investigated further the involvement of hsa-mir222-5p as one of the key factors responsible for the switch of melanoma cells. We, therefore, knocked down hsa-mir222-5p, showing that IgR39 cells loose the capability to revert the expression of CSC's markers and acquire a less aggressive phenotype, as shown by measuring the migration capabilities of the cells in 2D wound healing assays and in 3D spheroids invasion assays.
Altogether, our results shed a new light on tumor plasticity and provide guidance for further targeted studies aimed at impairing it.

Overshoot in the Mean-Field Model
We developed a mean-field model in which CSCs are capable of unlimited proliferation [1,23], while CCs can only divide into two CCs for a limited number g of generations, after which they become senescent and die. Details and equations of the model are reported in the Materials and Methods section and in Appendix A.
The model was used to generate the dynamical behavior of the number S of CSCs, of the number N j of CCs with various age j, of the total number N = ∑ g i=0 N i of CCs and of the number m of inhibitor molecules. A first experimental observation that the model should reproduce is that the total number of cells S + N must be an increasing and unbounded function of time. A necessary and sufficient condition for this is that Once the condition is satisfied, we can compare with experiments the fraction c = S S+N of CSCs over the total number of cells. This quantity has been observed [1] to experience an overshoot ten days after sorting in a population initially devoid of CSCs. Furthermore, the amount of CSCs in the overshoot peak decreases with the fraction c(0) of CSCs remaining in the initial population, as shown in Figure 1b (red bars).
The set of parameters associated with an overshoot which are most similar to the experimental data, is given in Table A1, corresponding to dynamics of c(t) and to a dependence on c(0) described by the continuous lines in Figure 1b. The initial conditions for these simulations imply an initial number of total cells of 10 6 [1], with an initial fraction of sorted CSCs ranging from 0.001 (0.1%) to 0.01 (1%). CCs are assumed to be distributed among the g generations reflecting the steady-state distribution among generations associated with our model. We assume that is m(0) = 0 during the sorting procedure when the inhibiting molecule is washed out (although a small amount of m(0) does not change qualitatively the results). The overshoot results to be sizable if the sorting procedure leaves up to 1% of CSCs in the initial set of cells (see Figure 1c), in agreement with experiments [1]. Consequently, there is at least a set of parameters for which the model reproduces qualitatively the experimental data. Being a thorough exploration of the whole parameter space of the model computationally unfeasible, we started from those of Table A1 and studied the dynamics of the system varying one parameter at a time.  The value of the maximum c max concentration of CSCs during the overshoot and of its stationary value c * are displayed in Figure A1, as function of some of the parameters of the system. The size of the overshoot increases smoothly but considerably as a function of the switching rate σ and of the rate difference , and decreases as a function of the aging rate k d of CCs. The dependence on the other parameters is much weaker. These results suggest that the overshoot is overall a robust property of the model. The strong dependence on the switching rate σ indicates that phenotypic switching is an important ingredient to reproduce the experimental observations. Two important parameters seem to be k d and g, that control together the total amount of CCs. For large values of k d or g, the overall number of CC is large and the overshoot effect disappears. A simple qualitative explanation for this is that if we increase the values of k d or g, then the total number of CC cells increases exponentially, while the number of CSCs increases only as there are more CC cells able to switch. Therefore we expect the increase in those values to reduce the amount of the overshoot.

Inhibition of Phenotypic Switching
Since the phenotypic switching is inhibited by a molecular species m 0 produced by CSCs, it should be possible to inhibit phenotypic switching after sorting out CSCs by growing the cells in a medium containing the inhibiting molecule. To understand this idea from a quantitative point of view, we performed simulations according to Equation (3) and compared the total number of cells at fixed time T when the sorted cells are treated at time t = 0 with a solution containing m 0 inhibiting molecules, with the total number of cells in the untreated case (i.e., m 0 = 0). We defined the quantity to quantify the effect of the inhibiting molecule at time T. It should be noted that the amount of inhibiting molecule varies with time from its initial value, because, besides that provided at time zero, it is produced by CSCs and undergoes degradation.
In Figure A2a, we show the variation of R T,m 0 as a function of the time T after sorting. The figure shows that a noticeable effect can be seen already a few days after sorting, although a strong effect could be seen after 20 days, when the number of cancer cells is reduced by about 50% for the highest initial concentrations of the inhibiting molecule. In Figure A2b, we also display the value of R 20days,m 0 as a function of m 0 , for different choices of c 0 , evaluating the number of cells after 20 days. Independently of the remaining concentration c 0 of CSCs, there is a threshold of ∼10 5 molecules, beyond which one can observe an inhibitory effect. Assuming that the plate used for the experiment has a volume of ∼10 mL, this threshold corresponds to a concentration of 10 −17 mol, meaning that even a tiny amount of molecule has a detectable effect on the (exponential) growth of the cells. The effect of the inhibitor saturates to a c 0 -dependent value when its copy number reaches ∼10 12 , corresponding to ∼1 nmol.

Simulations of Tumor Growth in Two Dimensions
The rate Equations (3) average in space the concentration of cells and of the inhibitor factor. To investigate whether specific spatial pattern appear during tumor growth, we studied an extension of the model in which the spatial degrees of freedom appear explicitly.
We assume that cells sit on the vertexes of a square lattice, whose elementary length is then ∼10 µm and whose linear size is typically 10 3 ∼1 cm. Each site can host at most one cell, either normal or CSC. It undergoes the same processes (symmetric or asymmetric division, phenotypic switching, cellular aging and death, production of the switching-inhibitor molecule) as described in Figure 1a for the rate-equation model. Also the numerical values assigned to the rates of the different processes are the same. The cellular dynamics is simulated with a Gillespie algorithm [24]. Upon cell division, one of the new cells remains in the same site, while the other occupies a neighboring site. If all the neighboring sites are occupied, all cells in a random direction are shifted to make room for the new one. Neighbors are defined along the horizontal, vertical and diagonal directions, according to the Moore scheme, which minimizes artifacts due to the square lattice [25]. A further parameter which was not defined in the mean-field model is the diffusion coefficient of the inhibiting factor. Assuming that it is a nm-sized miRNA complex, Stokes' law would suggest a diffusion coefficient D∼100 µm 2 /s, which would decrease to D∼0.1 µm 2 /s in case it is encapsulated into a vesicle.
Simulating the diffusion of the inhibiting factor within the Gillespie algorithm is computationally unfeasible, because each movement of every single molecule takes place on a time scale which is much smaller than that associated with the cellular processes, and thus essentially all the computational time would be used to simulate microscopic diffusion processes. Consequently, diffusion is simulated in parallel to the Gillespie algorithm by solving a discretized version of the diffusion equations as discussed in the Materials and Methods section. Simulation of the dynamics of this model starting from a set of cells purified from CSCs display an overshoot in CSCs concentration, similar to what was observed in the rate equations and in the experiments.
We then focused our analysis on the growth of a tumoral mass starting from a single CSC, as displayed in Figure 2a. The overall growth of the number of CCs and CSCs becomes exponential after ∼5 days with a rate that depends mildly on the diffusion coefficient of the molecule (from 0.4 days −1 at low D to 0.5 days −1 at large D) and is indistinguishable from a model without phenotipic switching. The fraction S/(S + N) of CSCs decreases to a stationary value of few percent with a weak dependence on D. Unlike the cellular growth rate and the fraction of CSCs, the spatial organization of CCs and CSCs depends markedly on the properties of the diffusing molecule. Figure 2a shows that for higher values of the diffusion constant, the clusters of CCs are more dense. The most notable feature of is the accumulation of CSCs at the border of the tumor. This effect increases with the diffusion coefficient D of the molecule and is absent if no switching applies. This effect can be better quantified from the radial distribution of the fraction of CSCs displayed in Figure 2b. When D is large, the fraction of CSCs at the border of the tumor displays a higher peak than at smaller values of D, while no enrichment at all is observed in the absence of switching.
From the snapshots displayed in Figure 2a it is apparent that for small and intermediate values of D, CSCs tend to cluster together, especially at the border. To get a more quantitative insight in this, we performed a clustering analysis of CSC based on their spatial density [26]. In Figure 2c we show the multiplicity of clusters found and the distribution of their density. Without switching, the system displays a small number of high-density clusters. A similar behavior is observed with switching at low diffusion constant. Increasing D, we observe many low-density clusters, suggesting that CSCs decrease their tendency to partition in clusters.
Summing up, clustering seems to appear as a consequence of CSCs duplication, while switching, which is more effective at large D, tends to make the distribution of CSCs more uniform. This is not unexpected, since cell duplication is a local phenomenon, while switching is random in space, and can only be modulated in space by the diffusing factor. The high diffusivity of the factor makes, however, its distribution to be quite uniform, suppressing most spatial patterns.

Impact of hsa-mir222-5p on Phenotopic Switching
According to the model when CSCs are sorted out by flow cytometry, CCs are not exposed anymore to the switch inhibiting factor and therefore start to massively switch in agreement with our previous observations [1]. In order to better understand the mechanism underlying the model hypothesis, we collected the conditioned media (CM) of IgR39 cells and analyzed the level of expression of miRNAs. As shown in Figure A3, there is a good correlation between the expression of miRNAs inside the cell and the amount of miRNAs found in the conditioned media. Since hsa-mir222-5p was strongly involved in phenotypic switching [1], we knocked it down in IgR39 cells by pEGFP-Sponge-mir222-5p (sh-mir222, see Materials and Methods). As shown in Figure A4, the level of expression of hsa-mir222-5p was significantly lower with respect to cells infected with a scramble miRNA lentivirus. Moreover under the same conditions, a phenotypic marker of CSCs, CXCR6, was expressed at very low level, see Figure A5. Finally, we checked if the knock down of hsa-mir222-5p had an impact on other miRNAs. As shown in Figure A6, no significant changes in expression were detected in other miRNAs (see Materials and Methods for more details).

Impact of hsa-mir222-5p on Cell Migration
Epithelial-mesenchymal transition (EMT) is a critical feature that helps tumor cell spreading by acquiring a higher motility [27]. Herein, we checked if the lack of expression of hsa-mir222-5p had an impact on cell migration using two different approaches: (a) in a 2D wound healing assay and (b) in 3D spheroids quantifying the invasion of the cells grown as spheroids into a collagen network [28]. In both 2D or 3D models we found a decrease in motility of the cells without hsa-mir222-5p, confirming the role of this miRNA in mesenchymal phenotype [1]. In particular, Figure 3a shows two snapshots of the local velocities and vorticities obtained by the PIV method where arrows represent velocities and color indicates vorticity. The snapshots show that upon hsa-mir222-5p knockdown, velocities tend to decrease. This is confirmed by computing the distributions of velocities that we report in Figure 3b.
Next, we carried out experiments growing the cells in spheroids embedded in collagen networks according to the protocol discussed in the Materials and Methods [28]. Time lapse observations are summarized in Figure 3c for scramble and sh-mir222 knockdown cells (see also Videos S1-S2 and Appendix C). We can see that while scramble cells spread into the matrix, sh-mir222 knockdown cells remain confined into the spheroid. While cells do not spread, we observe motion of cells within the spheroid (see Video S2). This implies that while motility is still present, invasion capabilities are strongly impaired. We quantified these observations by measuring the typical spread of the spheroid as ( (R − R c ) 2 ) 1/2 , where R c is the center of mass of the spheroid. Figure 3d summarizes the results.

Discussion
Melanoma is a highly resistant tumor with poor patient survival due to its intrinsically high metastatic capability [27]. The appearance of resistance after chemotherapy or the presence of intrinsic resistance leads to great difficulties in devising an effective and durable therapy. The ability of tumor cells to change their status using epigenetic mechanisms independent of the environment, such as the tumor niche, has been shown to play a critical role in the acquisition of a resistant phenotype in response to specific drugs. In light of these findings, the ability of tumor cells to change their phenotype, becoming resistant, can happen at the level of the tumor cells and/or at the level of the tumor niche. Interfering with this mechanism appears to be a promising and crucial issue to fight metastasis [27].
Our group recently demonstrated that human melanoma cells can change their phenotype, expressing EMT markers dynamically and activating Wnt-β-catenin pathway, thanks to a complex network of miRNAs which includes hsa-mir-222-5p [1]. The most important consequence of those findings is that tumors can control the number of CSCs in the bulk. In fact, these factors drive the switch of CCs into CSCs under specific environmental conditions [1,27]. It is clear that one interesting strategy could be to interfere with this process, finding a way to stop phenotypic switching to CSCs. This was the main goal of the present study. To better understand the regulation of the number of CSCs in the tumor cell population, we resorted to mathematical models. We first considered a mean-field model subdividing cancer cells into different compartments corresponding to CSCs and CCs at different generations. While CSCs can self-renew or give rise to CCs and CCs divide only for a finite number of generations, we considered the possibility for CCs to revert into the CSC state. This possibility is modulated by an inhibiting factor released by CSCs. With these assumptions, the model is able to faithfully reproduce the phenotypic switching observed experimentally [1].
We have then investigated a putative candidate for the inhibiting factor: hsa-mir-222-5p. This miRNA is strongly involved in phenotypic switching, being highly expressed before the overshoot of CSCs [1]. Here, we found that this factor is also present in the conditioned media (CM), in agreement with recent observations in melanoma [29]. Knockdown of hsa-mir-222-5p has important consequences for tumor aggressiveness: (i) the cells no longer express CSC markers and (ii) their migration capability is impaired.
We can reconcile our model with our observations considering a feedback loop driven by hsa-mir-222 (see Figure 4). In the steady-state, hsa-mir-222 is expressed by the cells and released outside, acting on paracrine cells. We reported in Ref. [1] that in unsorted melanoma cells, the level of expression of LEF1 shows a two fold increase upon silencing of hsa-mir-222 ( Figure 4). LEF1 is a key downstream target of the Wnt signaling pathway. The canonical Wnt signaling is activated when the Wnt ligand binds the Frizzled and LRP5/6 receptors. This leads to the stabilization of β-catenin which is transferred from the cytoplasm to the nucleus where it activates LEF1 (see Figure 4). It is known that MITF, which is directly induced by LEF1, acts as an activator of hsa-mir-222 [30]. Hence, activation of the Wnt pathway leads to increased hsa-mir-222 expression. On the other hand, it has been reported that hsa-mir-222 can inhibit the translation of Frizzled [31,32], thus suppressing the activity of the pathway (Figure 4). Accordingly, as already shown in [1] and reported in Figure 4, CCs three days after sorting, at the onset of phenotypic switching, display a level of LEF1 similar to the one in unsorted cells. When hsa-mir-222 was silenced, LEF1 decreased slightly. Taken together, these results suggest that the switch inhibitor is hsa-mir-222 itself within the feedback loop illustrated in Figure 4.
In addition to mean-field reaction equations, we have also studied a two-dimensional variant of the model where we considered the interplay between tumor growth, phenotypic switching and the diffusion of the switch inhibitor. Simulations of the model showed when phenotypic switching is present, CSCs cluster at the boundary of the tumor mass in a diffusion-dependent manner. This result is interesting because it provides a natural mechanism that could help the dissemination of CSCs into secondary tumors [27]. Notice that in a previous spatial model of cancer growth with CSCs, the localization of CSCs at the border was imposed by reducing the motility of the other cells [33,34]. The enrichment of CSCs at the border is apparent in the simulations both at large and small values of the diffusion coefficient D, although it is more marked in the former case. However, the reason for it seems to be slightly different in the two cases. At large D, the concentration of molecule is diluted out fast to the surrounding of the cellular system, displaying a radial profile that decreases moving away from the center and leaving the border with a concentration low enough to allow phenotypic switching. On the other hand, at low values of D, a transient initial enrichment propagates in time and in space along the radial direction, and is maintained at all later times. The diffusion constant D might play an important role also for possible therapeutic strategies based on mir222, since the size of carrier would determine the value of D and consequently the dynamics of phenotypic switching.
In conclusions, our paper elucidates the molecular mechanism of phenotypic plasticity in melanoma cells and provides guidance for therapeutic interventions targeting the switching into CSCs. To this end, it would be interesting to corroborate our results in vivo following the diffusion of fluorescent-tagged sh-mir222 and tracking the location of CXCR6 positive cells by immunohistochemistry.

The Computational Model
We construct a mean-field model in which CSCs are capable of unlimited proliferation [1,23] and can undergo different kinds of cellular division, into two CSCs, into two CCs, or into one CSC and one CC. On the other hand, CCs can only divide into two CCs for a limited number g of generations, after which they become senescent and die [23]. We also assume that CCs can switch stochastically to CSC at constant rate [35]. The new hypothesis of the current model with respect to previous models is that CSCs produce a molecular species that inhibits switching of CCs into CSCs according to a Michaelis-Menten molecular mechanism. We will discuss the experimental evidence for this in the second part of the paper.
We first introduce a delay in the action of CCs, representing the cellular replication and switch according to experimental data [1]. Under well-mixed microenvironment hypothesis, the model can be expressed by a system of delay rate equations in the concentration S of CSCs, the number N j of CCs at generation j (with 0 ≤ j ≤ g) and the number m of the inhibiting species. The overall mechanism is sketched in Figure 1a and its dynamics can be described by the following system of equations with initial conditions S(0) = S 0 , m(0) = m 0 and the N j (t) = N 0 for t ≤ 0 are set to the same initial values for sake of simplicity. The rate equations are solved numerically using a custom made C code. The numerical method used to solve the system of differential equations is a variable-coefficient linear multistep Adams method in Nordsieck form (from the GNU Scientific Library), which employs a variable time-step.
Some of the numerical values that define Equation (3) are known, at least as order of magnitude. The division rate of melanoma cells add the variable, as in × 0.5-0.6/day. Typical initial populations of CCs is typically N 0 ∼10 6 and that of CSCs is approximately one hundredth of it [1,23,36].
For the two-dimensional model we solve a discretized version of the diffusion equations where m(r, t) is the number of molecules at position r of the lattice at time t and ∆t = 3 × 10 −4 days is the time step used for numerical differentiation, following a five-point stencil finite-difference method.

Conditioned Media
Cells were growth to 90% confluence and maintained for 16 h in the standard growth medium condition without 10% FBS. The media were collected and after a brief centrifugation to discard dead cell and debris, the supernatant was stored at −80 • C for no more than 1 month until use.

FACS Analysis
Cells were analyzed for phycoerythrin (PE) anti-human CXCR6 (cod. MAB699, R&D System, Minneapolis, MN, USA). Non-specific mouse IgG is used as isotype control (Isotype). Non-specific mouse IgG2B (R&D, Minneapolis, cod.IC004IP) is used as isotype control (Isotype). For each flow cytometry evaluation, a minimum of 10 6 cells were stained and at least 50,000 events were collected and analyzed. Flow cytometry and analysis was performed using a Gallios flow cytometer (Beckman Coulter Indianapolis, IN, USA).

miRNA Analysis by miRCURY LNA Universal RT miRNA PCR
Total RNA was extracted using the standard RNA extraction method with TRIzol TM (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA). RNA concentration in each sample was assayed with a ND-1000 spectrophotometer (NanoDrop, Thermo Fisher Scientific, Waltham, MA, USA). Reverse transcription reactions were performed using the miRCURY LNA RT Kit (cod. 339340, Exiqon, Qiagen, Germany). Reverse transcription thermocycling parameters were as follows: 42 • C for 60 min and 95 • C for 5 min. Real time q-RT-PCR analysis was performed using ViiA 7 Real-Time PCR System (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA). The PCR reaction included 1 ng of template cDNA, 1 µL of LNA PCR primers mix, 1 µL of RNAse-free water and 5 µL of miRCURY SYBR Green Master Mix (cod. 339345, Exiqon, Qiagen, Germany) in a total volume of 10 µL. Cycling conditions were as follows: 95 • C enzyme activation for 10 min, followed by 50 cycles of amplification: 10 at 95 • C for denaturing, 1 min at 60 • C for annealing/elongation. Melting curve analysis was performed between 60 • C and 95 • C at a ramp rate of 0.11 • C/s.

Isolation and Detection of miRNAs
Total RNA was extracted using the standard RNA extraction method with TRIzol TM (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA). RNA concentration in each samples was assayed with a ND-1000 spectrophotometer (NanoDrop, Thermo Fisher Scientific, Waltham, MA, USA) and its quality assessed with an Agilent 4200 TapeStation or Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Next generation sequencing experiments on infected cell lines was performed by Genomix4Life (Baronissi, Salerno, Italy). Indexed libraries were prepared from 1 µg/ea purified RNA with TruSeq SmallRNA Library Prep Kit (RS-200-0012, Illumina, San Diego, CA, USA) according to the manufacturer's instructions. Libraries were quantified using the Agilent 4200 TapeStation (Agilent Technologies, Santa Clara, CA, USA) and pooled such that each index-tagged sample was present in equimolar amounts, with a final concentration of 2 nM for pooled samples. The pooled samples were subject to cluster generation and sequencing using an Illumina NextSeq 500 System (Illumina, San Diego, CA, USA) in a 1 × 75 single end format. The raw sequence files generated underwent quality control analysis using FastQC [38].

miRNA Differential Expression Analysis
Relative expression values of miRNAs are computed taking the average of the scramble as reference values. We filtered the initial pool of 2656 miRNAs as follows: First, we discard miRNAs with absolute expression below 10 3 in any of the samples. Then, we keep only those miRNAs whose expression displays changes of 10% with respect to the Scramble samples. Finally, we keep only those miRNAs that display a change of at least 25% with respect to Scramble in at least one of the Sponge samples.

miRNAs Expression in Conditioned Media
MiRNAs were extracted from the conditioned media collected as described in the materials and method section by miRCURY RNA Isolation kit-Biofluids (cod. 300112, Exiqon, Qiagen, Germany). According to the manifacturer's instructions, the samples were treated with DNase and eluted with 20 µL nuclease-free water. The quality of RNA extracted was carried out by capillary electrophoresis using Bioanalyzer (RNA 6000 Pico Kit, Agilent Technologies, Santa Clara, CA, USA). The libraries were obtained using the kit TruSeq Small RNA Library Prep (RS-200-0012, Illumina, San Diego, CA, USA) and validated by Bioanalyzer using Agilent High Sensitivity DNA kit (cod. 5067-4626, Agilent Technologies, Santa Clara, CA, USA). A pick of about 140 bp was expected corresponding to miRNA, since during the preparation of the libraries we added two 60 bp adapters to each end of the existing 22 bp miRNAs sequence. To enrich the libraries in the fraction of the miRNAs only, we carried out an 6% polyacrilamide electrophoresis (gel Novex TM TBE, Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA) and we cut the band between the two markers of 145 and 160 bp. We then validated the libraries and quantified by Bioanalyzer using Agilent High sensitive DNA (cod. 5067-4626, Agilent Technologies, Santa Clara, CA, USA). The library was then sequenced by Genomix4Life (Baronissi, Salerno, Italy) using Illumina NextSeq 500 System (Illumina, San Diego, CA, USA).

Migration Assay
A wound was introduced in the central area of the confluent cell sheet by using a pipette tip and the cell migration was followed in time-lapse using a Leica DM IL LED inverted microscope (Leica, Germany) at 37 • C and 5% CO 2 and 95% humidity (UNO stage top incubator, Okolab, Italy) every 15 min until the complete closure of the scratch with a 10× objective.

PIV Analysis
Velocity field was estimated using PIVlab app for Matlab [42]. The method is based on the comparison of the intensity fields of two consequent images. The difference in the intensity is converted into velocity field measured in px/frames and then converted to µm/h [43].

Spheroids Culture, Collagen Invasion Assay and Image Analysis
Multicellular spheroids were obtained from sub-confluent cells using the hanging-drop technique (2% methylcellulose, cod. M6385, Sigma Aldrich, Merck, Darmstadt, Germany; 500 cells/50 µL drop; 7 days spheroid assembly) according to Ref. [44]. Spheroids were embedded in non-pepsinzed rat-tail collagen type I solution (2.5 mg/mL; BD Biosciences, Franklin Lakes, NJ, USA) prior to collagen polymerization (37 • C for 5-10 min). The collagen matrix invasion process was monitored by bright-field time-lapse microscopy using a LEICA DMI8 microscope (Leica, Germany) acquiring one frame per hour for 72 h. The image analysis is done using Matlab Image Processing Toolbox. At first, the image is cut to focus on the area of interest. Then it is thresholded and a black-white image is created using built-in Matlab function "edge". After it is dilated with "imdilate" and the holes are filled with "imfill". Small noise is removed by "bwareaopen". Then the image is eroded and filled by "imerode" and "imfill" functions. In the resulting images the center of the white (detected) area is found. Then, we compute the average distance from the center to all other detected points. This value corresponds to the spread at a given time step.

Statistical Analysis
For miRNA expression, p-values were obtained from EdgeR differential analysis [45]. For all the other cases, the two-tail unpaired t-test was used to assess statistical significance.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

CSC
Cancer stem cell CC Cancer cell SC Stem cell sh-miR222 Sponge miR222

Is the Mean-Field Model a Minimal Model?
A relevant question about the mean-field model is whether all the terms in Equation (3) are really necessary to describe the system. If the time delay τ is not introduced, the overshoot peak always occurs within the few hours after sorting, while it is observed that some days are needed for the increase of c(t) [1]. If one disregards the aging of CC, describing them just in terms of their total number N, transforming Equations (3) into then it is possible to write a nearly closed equation for the concentration c(t) of CSC. Setting ∆ = k h k h +m h , this reads Since for large time and > 0, then m −→ ∞, and so ∆ = k h k h +m h −→ 0, giving the stationary value For c * to be of the order of the observed order of magnitude ∼10 −2 , Equation (A2) requires that 3k 2N + k 1N1S k 2S + k d . This is unlikely to be the case, because the experimentally observed division rates k d , k 2S , k 2N , k 1S1N are all of the order of 0.1 days −1 [23] , and for sure k 2N < k 2S , because we are in the case > 0.
Finally, if phenotypic switching occurs in an uncontrolled way, no overshoot is possible. In fact, if one removes the dependence on m in the second of Equation (A1), this is equivalent to set ∆ = 1 in Equation (A2). Setting f (c) = −( + λ + k d c 2 + (k d + )c), the stationary points of Equation (A2) can be found from f (c * ) = 0 and are These are always real because + λ + k d > 0, and have opposite sign (let us call, e.g., c * 2 < 0 and c * 1 > 0). The solution c * 2 is not admissible because it cannot describe a fraction of cells. The stability of c * 1 is determined by the condition which is satisfied for every initial concentration c 0 in the range [0, 1]. One can exclude the existence of an overshoot showing that c(t) is monotonous as it approaches c * 1 for every t and every c 0 . Since −( + λ + k d < 0, then f (c) > 0 if and only if c * 2 < c < c * 1 , that for the case of interest becomes 0 < c(t) < c * 1 . For c(t) to have an overshoot, starting from a small c 0 , c(t) must grow in a finite time to a peak strictly above the attractive equilibrium value c * 1 . But this is not possible since c (t) = f (c) is strictly negative in the region c * 1 < c(t) < 1. Thus, in the case of undamped (non-inhibited) phenotypic switching at constant rate σ the phenomenon of overshoot of the concentration c is not possible.    Figure A4. Efficiency of hsa-miRNA-222-5p knockdown in IgR39 cells. MicroRNAs silencing efficiency was evaluated by measuring hsa-miRNA-222-5p in IgR39 sponged cells (sh-mi222 5-p) according to the Materials and Methods section by relative expression and standard errors over three replica with qRT-PCR (a) (** p < 0.005) and microRNAs abundance (b). Significance has been assessed using EdgeR (* p-value < 0.05). Cells sponged with a scramble miRNA sequence (Scramble, see the Materials and Methods section) was used as a control. sh-mir222 Scramble Figure A6. Expression of miRNAs 10b-5p, 15b-5p, 221-3p and 378a-3p for WT (a) and Sponge hsa-mir222-5p (b) IgR39 cells, compared to the Scramble condition, see Materials and Methods section for details. All comparisons lead to non-significant p-values when using a two-sample t-test. Table A1. Parameters of the model and corresponding molecular/cellular processes.  Table A2. Sense and antisense primers for seqeuncing of pEGFP-Sponge vector used for verify the presence of 8 inhibitor sequences for hsa-mir222-5p according to the manufacturer's instructions.