Drugs modulating stochastic gene expression affect the erythroid differentiation process

Summary

To better understand the mechanisms behind cells decision-making to differentiate, we assessed the influence of stochastic gene expression (SGE) modulation on the erythroid differentiation process. It has been suggested that stochastic gene expression has a role in cell fate decision-making which is revealed by single-cell analyses but studies dedicated to demonstrate the consistency of this link are still lacking. Recent observations showed that SGE significantly increased during differentiation and a few showed that an increase of the level of SGE is accompanied by an increase in the differentiation process. However, a consistent relation in both increasing and decreasing directions has never been shown in the same cellular system. Such demonstration would require to be able to experimentally manipulate simultaneously the level of SGE and cell differentiation in order to observe if cell behavior matches with the current theory. We identified three drugs that modulate SGE in primary erythroid progenitor cells. Both Artemisinin and Indomethacin decreased SGE and reduced the amount of differentiated cells. On the contrary, a third component called MB-3 simultaneously increased the level of SGE and the amount of differentiated cells. We then used a dynamical modelling approach which confirmed that differentiation rates were indeed affected by the drug treatment. Using single-cell analysis and modeling tools, we provide experimental evidence that, in a physiologically relevant cellular system, SGE is linked to differentiation.

Keywords:

Gene expression – Cell differentiation – Drug therapy – Drug discovery – Entropy – Drug interactions – Artemisinin – Kullback Leibler divergence

Introduction

Cell-to-cell variability is intrinsic to all living forms, from prokaryotes [1, 2] to eukaryotes [3]. Such variability originates from many sources, but arguably stochastic gene expression (SGE) is an important driving force in the generation of cell-to-cell variability among genetically identical cells [4], although additional regulation layers do exist [5]. Classically, SGE is separated into intrinsic and extrinsic sources [6–10] even if in many cases distinguishing between the two is difficult.

The very existence of SGE led to the concept of a probabilistic mapping between inputs (environment) and outputs (cell decisions) [11]. It is therefore clear that SGE has to be precisely tuned so as to tailor the biological process in which it is involved [12, 13].

As an inevitable consequence of the molecular nature of gene expression process, it is clear that mechanisms dedicated to its control must exist. For example, it has been discussed that during the development, to ensure the reliable reproduction of some particular patterns, cell-signaling devices (e.g. the Wnt pathway) act as a noise-filter [14].

By contrast, numerous arguments suggest that SGE variation plays an important role in a wide range of biological processes ranging from bet hedging [15] to the fractional killing of cancer cells [16]. The involvement of SGE in decision-making has been shown in viruses [17–19], in prokaryotes [20–22] and, more recently, in metazoan cells [22–24]. Evidences of its role in the differentiation ability began to emerge from literature through different points of view.

First, the theoretical arguments [25–27], based on evidence of cell-to-cell heterogeneity but without any experimental demonstration, represent differentiation seen from a dynamical point of view [28]: undifferentiated cells exhibit a stable gene expression pattern corresponding to a pre-existing self-renewal state. Once differentiation is activated through external factors, cells could increase their SGE, explore a broader region of their state space and exhibit a mixed gene expression identity between the original and the destination states. Such an exploratory behaviour would increase the probability for cells to attain the space region where they stabilize their new gene expression pattern by reaching another pre-existing stable state: the differentiation state [26, 27]. In other words, an increase of SGE would lead to an improvement of the differentiation dynamic [26]. In addition, it has been described that, in a transition state, the original state may be removed as the gene expression values change, making the process irreversible [27, 29, 30]. Collectively, it has been largely discussed that stochastic fluctuations, associated with extrinsic inputs, appear to be the main means by which cells can change their state [22, 27, 30]. However, this theoretical point of view needs to be supported by biological observations.

Secondly, observation of the role of SGE during differentiation emerged, some years later, through experimental measures of the amount of SGE during differentiation processes [31–34]. We recently described a surge in cell-to-cell variability that precedes an irreversible differentiation of normal primary chicken erythroid progenitors called T2EC [35], that is fully compatible with such a view [31]. Interestingly, these results have been confirmed in various settings, ranging from the differentiation of murine lymphohematopoietic progenitors [33] to the differentiation of murine embryonic stem cells [32, 34].

Finally, experimental demonstrations such as studies of artificial modulation of the amount of SGE [13, 36] started to emerge. These last types of evidences aim at reinforcing the consistency of the relation between SGE and differentiation but a clear demonstration and characterization of this link is still lacking. Here, we pave the way toward this demonstration, adding a new complementary study to the legacy.

To do so, it is necessary to show that pharmacological modulators such as drugs [12, 13, 36] would on one hand modify SGE and, on the other hand modulate the differentiation process. It has recently been described that such drugs, identified using a large screening approach, were able to modulate the noise affecting a LTR-driven reporter gene in human T-lymphocytes [37]. In addition, drugs that directly inhibit promoter nucleosome remodelling were also shown to provide fine-tuning of SGE [38].

In order to demonstrate the general aspect of the relation between SGE and differentiation, we decided to explore the extent to which some of those drugs, that are able to reduce (Artemisinin and Indomethacin) or to increase (MB-3, [7]) the level of SGE, could alter differentiation.

Here we show that the three selected drugs modify, significantly and simultaneously, the level of SGE and the rate of cell differentiation. We therefore provide a clear evidence that, in a physiologically relevant cellular system, a pharmacological modulation of SGE is positively and consistently accompanied by a modification of differentiation, as suggested by existing points of view.

Materials and methods

Cell culture and treatment

T2EC were extracted from the bone marrow of 19 days-old SPAFAS white leghorn chickens embryos (INRA, Tours, France). These cells were maintained in a medium called LM1. It is composed of α-MEM medium supplemented with 10% Foetal bovine serum (FBS), 1 mM HEPES, 100 nM β-mercaptoethanol, 100 U/mL penicillin and streptomycin, 5 ng/mL TGF-α, 1 ng/mL TGF-β and 1 mM dexamethasone as previously described [35]. T2EC were induced to differentiate by removing the LM1 medium and placing cells into the DM17 medium (α-MEM, 10% foetal bovine serum (FBS), 1 mM Hepes, 100 nM β-mercaptoethanol, 100 U/mL penicillin and streptomycin, 10 ng/mL insulin and 5% anemic chicken serum (ACS)). Differentiation kinetics were obtained by collecting cells at different times after the induction in differentiation. For Indomethacin and Artemisinin, cells in the self-renewing medium are treated at respectively 25 μM and 1 μM 48h before switching into a differentiated medium in order to optimize their effects. For MB-3 ((2R,3S)-rel-4-Methylene-5-oxo-2-propyltetrahydrofuran-3-carboxylic acid), cells are treated at 10 μM just after inducing the differentiation. For each drug, a control treatment (0.1% DMSO) was added following the same conditions.

Counting of cell viability and cell differentiation

Cell population growth was evaluated by counting living cells using a Malassez cell and Trypan blue staining (SIGMA). This method was also used to assess the toxicity of the drugs in T2EC (S6 Fig). Cell population differentiation was evaluated by counting differentiated cells using a counting cell and Benzidin (SIGMA) staining which stains haemoglobin in blue.

Dynamical model for erythroid differentiation

Every detail regarding the design, selection, calibration or identifiability analysis of our dynamic model can be found in its original paper [39]. All data and pieces of code used for the curent study are available in a public github repository https://github.com/rduchesn/Drugs-modulating-stochastic-gene-expression-affect-the-erythroid-differentiation-process. Herein, we give only the definitions of the useful concepts and methods in this study.

Model definition

Dynamic model The ODE governing the time-evolution of the cell populations in each compartment of the model are given in Eq 1:

It is characterized by the set (ρ_S, δ_SC, ρ_C, δ_CB, ρ_B) of five parameters, where ρ_i is the net proliferation rate of compartment i (positive or negative), and δ_ij is the differentiation rate of cell type i into cell type j, which is positive.

The solution to this linear ODE model is given in the supplementary materials of the original paper [39]. We will note f_i(t_j, y₀, θ) the prediction of this model for the i^th observable of the experiment, on the j^th timepoint t_j, with the initial condition y₀ and parameters θ.

Error model In order to properly define the statistical likelihood of our model, we introduce a gaussian model for the distribution of the residuals of the dynamic model:

Estimation in the control case

Likelihood From Eq 2, the likelihood of the model naturally follows, and we can estimate the best-fit parameter values of the model by minimizing the negative logarithm of the likelihood:

Algorithmic details We minimized −2 log(L) using the Truncated Newton’s algorithm [40, 41] implemented in the python package for scientific computing scipy [42]. Convergence to the global minimum was assured by a random sampling of the initial guesses for parameter values.

Estimation under treatment

Parameter variations The model has seven parameters (five dynamic parameters: ρ_S, δ_SC, ρ_C, δ_CB and ρ_B; and two error parameters: b₁ and b₂), of which six are estimated from the data. Under a given treatment, we consider that each estimated parameter could either be estimated from the data, or set equal to its control value. For each treatment, this defines 2⁶ = 64 different models, with a varying number of additional parameters to estimate.

Model selection We estimated the parameter values of these 64 models for each treatment, and selected the best models by computing their corrected Akaike’s Information Criterion [43]:

Single cell high-throughput RTqPCR

Every experiment related to high-throughput microfluidic-based RT-qPCR was performed according to Fluidigm’s protocol (PN 68000088 K1, p.157-172) and recommendations. All the following steps from single-cell isolation to high throughput RTqPCR of each cells are described in [31].

Entropy

We estimated the Shannon entropy of each gene j at each timepoint t as follows: we computed basic histograms of the genes with N = Nc /2 bins, where Nc is fixed for all tests, which provided the probabilities p j , k t of each class k. Finally, the entropies were defined by

Results

Drugs affect noise

In order to characterize the relation between SGE and differentiation, we first make sure to be able to change the amount of SGE in T2EC using three drug treatments: Artemisinin, Indomethacin and MB-3.

Artemisinin and Indomethacin are known to modify SGE of the HIV LTR promoter in human T-lymphocytes [37]. MB-3, a chromatin modifier, is known to modify stochastic gene expression in yeast [7] and in murine embryonic stem cells [13]. At first, we wanted to confirm that these drugs do indeed modify SGE in our cellular system and to determine the mechanisms associated with this effect.

We treated T2EC with or without drugs and induced their erythroid differentiation. We then performed single-cell high-throughput RTqPCR on these cells at different time points after differentiation. We assessed a 92 gene panel, relevant for erythroid differentiation study, identical to those previously measured in untreated cells [31].

There are various ways of quantifying the amount of so-called “noise” in gene expression. Unfortunately, no consensus has emerged with authors advocating for the use of normalized variance (NV = σ²/μ²) and claiming that the use of Fano Factor (F = σ²/μ) might be misleading [44, 45], whereas others defend the exact opposite position [46]. The use of the coefficient of variation (CV = σ/μ) is also known to be limited [45, 47–49].

Both we [31] and others [32, 50, 51] have recently proposed the shannon entropy as an alternative measure, that is dedicated to quantify cell-to-cell variability. For our purpose and based on its mathematical definition, a deterministic pattern of expression exhibits low entropy whereas a high entropy indicates a more diverse expression pattern [31, 50]. We therefore analyzed our single-cell transcriptomic data using this metric.

We can observe in (Fig 1A) that the entropy was affected by all treatments. Under Indomethacin or Artemisinin treatment, entropy significantly decreased after 2 days of erythroid differentiation. This effect was more pronounced with Indomethacin. The opposite effect is observed with the MB-3-based treatment, for which entropy was significantly increased after 12h of differentiation (Fig 1A).

We then assessed whether the same genes exhibit entropy variations under the different drug treatments or not. To do so, we computed a correlation value between the entropy variations for each pair of drugs. If the same genes are affected by two drugs, then one would expect their entropy variations to be correlated. We observed a significant correlation only for the genes affected by Indomethacin and Artemisinin treatment. MB-3 treatment seemed to affect the variability of a different set of genes (Fig 1B).

The entropy variation could be achieved by modulating the global mean gene expression or the gene expression variance. Thus we finally wanted to test if our drug treatments affected entropy through the modulation of the mean gene expression value. If so, one might expect to see a correlation between the variation of entropy and the mean expression level under drug treatment.

Indeed for two drugs out of three, Artemisinin and MB-3, we observed a significant inverse correlation between mean and entropy (Fig 1C). Nevertheless the effect of Indomethacin on entropy was not related to an effect on mean gene expression.

Here we have found three drugs that modulate SGE in T2EC cells. Indomethacin and Artemisinin decreased it whereas MB-3 increased it. MB-3 involved a different set of genes than Artemisinin and Indomethacin (also shown in S1 Fig) and the effect of drugs was not strongly related to an effect on mean gene expression value. Entropy modulation is therefore the only common characteristic of our three drugs.

Even with a weak correlation, it seems that Artemisinin and Indomethacin could affect the same genes to change the level of SGE. Thus, we decided to get more insight on genes targeted by each drugs by selecting genes of interest in function of the influence of drugs on them (Table 1).

We observed that for some of the genes that we studied, the entropy was affected by both Indomethacin and Artemisinin. We therefore decided to focus our analysis more specifically on these genes.

From the Shannon entropy value of each gene for each condition, we calculated the relative value for each treatment compared to the control and reordered these results from the highest ratio to the weakest. For a given gene in a given treatment, if this relative entropy value is very different from one, it means that the entropy of that gene was greatly affected by the treatment. In such a way, genes are represented in function of their entropy ratio values from top to bottom as the most positively affected by drugs to the most negatively affected (Table 1). The gene expression distributions of the most affected genes (negatively and positively) for each drugs were represented in S2 Fig. The first quartile of the resulting distribution of relative entropy indicates genes for the entropy was the most positively affected by the treatment. Conversely, the third quartile indicates genes for which the entropy was the most negatively affected by the treatment.

We can observe that for Indomethacin and Artemisinin treatments, the 3 top genes that are the most positively affected by the drugs in variability are the same: ALAS1, TPP1 and MFSD2B. Whereas, for MB-3, the relative entropy of TPP1 and ALAS1 genes are negatively affected (under the third quartile). Such as Indomethacin and Artemisnin, MB-3 affected positively the entropy of MFSD2B. We can also find the gene MKNK2 as a common positively affected gene for Indomathacin and Artemisinin (all their relative entropy values are superior to the first quartile). Under MB-3 treatment, the relative entropy of MKNK2 does not seem to be affected as its value is between the first and the third quartile (Table 1). For the genes that are most negatively affected by drugs, the results are less clear. Under the third quartile, we can find MTFR1 and LDHA for both Artenisinin and Indomethacin. For MB-3, LDHA is negatively affected. HRAS1 is also negatively affected by Artemisinin and Indomethacin but its value is not under the third quartile for the Indomethacin treatment. HRAS1 is positively affected under the MB-3 treatment but for MTFR1, the relative entropy value is unchanged (Table 1). Regarding gene expression distributions, the most positively affected genes (in terms of entropy), show wider distribution in the treated condition than in the control (S2 Fig). The inverse is true with the most negatively affected genes. That supports the relevance of the indicator used, Shannon’s entropy. Moreover, all the distributions are heavy-tailed, which is expected from mRNA single cell distributions [52].

As we have shown previously, to modulate the level of SGE, drugs did not target a specific set of genes. This is clearer for MB-3 compared with the two others than between Indomethacin and Artemisinin. However, the increase of the noise level of ALAS1, TPP1 and MFSD2B seems to be a common effect of Artemisinin and Indomethacin.

In addition, to support these results, we decided to analyse in silico the connections between the three drugs’ targets.

To do so, we compared the different targets known in literature. For MB-3, the only target known is KAT2A protein [7, 13, 53]. Indomethacin targets both Cyclooxygenases (COX-1 and COX-2 also called PTGS for Prostaglandin-Endoperoxide Synthase) [54]. For Artemisinin, the task to find its targets is more complex because of the unspecificity of this drug [55]. In 2019, Heller and Roepe listed targets of Artemisinin-based drugs among three proteomic studies [56].

All together, known connections between proteins were represented using the STRING database (http://string.embl.de/) in (S1 Fig). Each edge between two proteins corresponds to a known association between those proteins. We can observe that KAT2A is alone in the connection network. Both PTGS-1 and PTGS-2 are highly correlated together but poorly correlated with the rest of the network. Each link refers to a co-mention between these two terms in a PubMed Abstract. None of them [57, 58] shows a direct interaction between these molecules. All these results suggest that there is no direct interaction between the drug targets reported in the literature. However, we have to keep in mind that we only compare the data reported in literature and that a potential interaction between drugs remains possible but not yet discovered.

We next used these drugs to test their effect on the erythroid differentiation process.

Drugs affect differentiation

In order to know if drugs modulating SGE also affect the differentiation process, we measured the percentage of differentiated cells in treated and untreated conditions during 96h of erythroid maturation (Fig 2).

A significant modulation in the percentage of differentiated cells was observed for all three drugs (Fig 2).

Indomethacin and Artemisinin decreased the percentage of mature cells from 48h of differentiation onward. MB-3 acted earlier: it significantly increased the percentage of differentiated cells by 24h before returning to somewhat below the control level.

Indomethacin and Artemisinin, two drugs that decreased SGE, reduced the percentage of differentiated cells. Inversely, MB-3 that increased SGE, enhanced the percentage of differentiated cells.

However, at this stage, we cannot conclude that a modification of the level of SGE by drugs is associated with a change of the differentiation process itself. Indeed, these effects might have several origins including modification in growth or death rates of our cells, which we cannot measure experimentally. To decipher between these effects, we decided to use a mathematical model describing the dynamics of the in vitro erythroid differentiation [39].

Cellular basis of drug effect

Our model describes the dynamics of three cell populations related to three different stages of differentiation [39]. The first one is the self renewing state (S) where differentiation has not started; the third one is the differentiated state (B) where cells have finished differentiating. The second one is the committed state (C), comprising intermediary cells that are committed to differentiation but not yet fully differentiated (Fig 3). The dynamics of these compartments follow a set of linear ODE. From the size of the cell population in the culture (S3 Fig), it seems reasonable to use a deterministic framework when modelling the growth and differentiation of the whole population.

Our dynamic model is characterized by a set of five parameters θ = (ρ_S, δ_SC, ρ_C, δ_CB, ρ_B):

• ρ_i is the proliferation rate of compartment i, involving the balance between cell proliferation and cell death. This value can be either positive (more proliferation than death) or negative (more death than proliferation).

• δ_ij is the differentiation rate of cell type i into cell type j, which is positive.

Considering that there remains no self-renewing cells after 2 days of T2EC differentiation (S4 Fig, [31]), δ_SC is a fixed parameter fully determined by ρ_S [39].

In order to get the best description of the drugs effects with the fewest parameters, we used the same approach as described in [39] and in the Methods section.

First, we estimate the parameters of the dynamic model in the control condition, using the data presented in Fig 2 (living cells and differentiated cells counts in the self-renewal and the differentiation media). We have already proven that our model identifiable, both theoretically and practically [39], using the profile likelihood approach [59]. It thus makes sense to compare the parameter values between the treated and untreated conditions.

For a given treatment, we consider that each parameter could either be equal to its untreated value, or to another value which should be estimated from the data (thus introducing a new parameter in the model). We test all the combinations of parameters that might vary under each treatment, and we select the best ones using Akaike’s weights, that are displayed on S5 Fig [43].

In the end, the parameter sets that we display in Fig 4 are those that reproduced well, and with the fewest additional parameters, the cellular kinetics during the in vitro differentiation (S3 Fig).

For all of those best models, the parameter values for each treatment are displayed in Fig 4. Under Indomethacin or MB-3 treatment, ρ_S (net growth rate of the immature cells) was not affected in all models and slightly decreased under Artemisinin treatment. Therefore, δ_SC was not affected by the treatments either, since its value is entirely determined by the value of ρ_S.

Concerning ρ_C, the net growth rate of the committed compartment, its values were reduced compared to the untreated condition for the majority of models under Indomethacin or Artemisinin treatment, whereas for MB-3 its value increased in all models.

A more variable change between drug effect was observed with parameter ρ_B, which describes the net growth rate of differentiated cells. Under Indomethacin treatment, some of the best models did not show a different value when compared to untreated condition whereas some models displayed a reduced parameter value. Under Artemisinin treatment this value was unchanged for four models among five and increased for the other one. With MB-3 treatment, ρ_B decreased in all models.

Finally, we found that the δ_CB parameter, representing the differentiation rate between committed compartment and mature cell compartment was affected by all three drugs: both Indomethacin and Artemisinin reduced this differentiation rate whereas MB-3 increased it in all best models.

These results demonstrate that all three drugs alter the differentiation process by modifying all dynamical parameters including the differentiation rate between committed and mature cells. It is clear that drugs that reduce SGE decrease the differentiation rate T2EC and inversely that the drugs increasing SGE accelerate cell differentiation, in line with our initial hypothesis.

Discussion

In this study, we highlighted for the first time the existence of a relation between the pharmacological modulation of stochastic gene expression and differentiation in both directions (increase or decrease) in the same cellular system. We first showed that three drugs, selected from literature [7, 37], increase or decrease the level of SGE in our cells. We therefore tested their effect on the differentiation ability of avian erythropoietic progenitors. We identified which differentiation parameters were affected by drugs using a dynamical model of the in vitro erythroid differentiation [39]. We demonstrated that the modulation of the differentiation process impacted the differentiation rate between the last two compartments. We therefore demonstrated that drugs modulating the amount of SGE simultaneously modify the differentiation process supporting all existing points of view [22, 25–27, 30] and reaching toward recent experimental evidences [13, 36].

Indomethacin, Artemisinin and MB-3 have clearly different functions. Artemisinin is an antimalarial drug used against a parasitic infection [53]. Indomethacin is an anti-inflammatory drug that affects the prostaglandin pathway [54]. These drugs were selected from another study [37] for their effect on the level of SGE on a HIV LTR promotor in human lymphocyte line. In this study, Artemisinin and Indomethacin increased the SGE amout of the LTR promotor. The opposite effect between T2EC and the LTR promotor system could have different origins:

• First, the LTR promotor is the only DNA region analyzed in the original study [37]. In our study, we analyzed the SGE variation of 92 genes previously selected to be relevant for avian erythropoiesis.

• Secondly, we analyzed the effect of Artemisinin and Indomethacin on SGE in T2EC as a sum of each effect on numerous genes. For some of them, the level of SGE was increased while for the others, the level of SGE was decreased.

• Finally, the cells used in the original study and ours are completely different. The original study used a line of human T-lymphocytes whereas we used a primary culture of avian erythroid progenitors.

MB-3 is an inhibitor of GCN5, a histone acetyl transferase (HAT) that activates global gene expression [60]. Even in such a seemingly well-defined case, it should nevertheless be remembered that a very complex relationship may lie between the biochemical action of a drug (HAT inhibition) and its biological effect on SGE [19].

Considering these different functions, it is hard to imagine that all these drugs have in common anything else than their ability to modulate SGE in T2EC. Nevertheless, it is important to note that it remains difficult to be certain that the effect of the drugs on differentiation is due to a direct effect of the drug on its know target or on some off-target effect. One method to resolve this issue would be to act on the pathway that the drug targets, without using the molecule (e.g. knockdown a drug’s target). For MB-3, a knockdown of KAT2A was performed in parallel of the use of MB-3 in mesendodermal differentiation (Moris et al., 2018). Both affected differentiation and the SGE in the same manner, demonstrating that the results indeed seem to be directly related to the biological inhibition of KAT2A and not to another independent pathway. Performing the same experiments for the two other drugs would clearly address this issue, although this might be more challenging for Artemisinin, for which many targets have been identified [56].

The question then arises of the mechanisms through which these different drugs modulate SGE. We first assessed whether these drugs affected the entropy of the same genes. For Indomethacin and Artemisinin, we showed that indeed the entropy of some of the same genes were affected but with a weak correlation. In contrast, MB-3 increased SGE through a different set of genes. This tends to indicate that the modulation of cell-to-cell variability per se, relatively independently of the gene function involved, is related to a modification of the differentiation process (see below).

Even if the correlation is weak, drugs reducing cell-to-cell gene expression variability seemed to affect the entropy of the same genes. Indeed, Artemisinin and Indomethacin both increased the cell-to-cell variability of 4 genes: ALAS1, TPP1, MFSD2B and MKNK2. ALAS1 gene encodes for a protein involved in the heme production in red blood cells [61]. TPP1 gene, previously named CLN2, encodes for a soluble lysosomal enzyme involved in metabolism [62]. MFSD2B is a lipid transporter released by erythrocytes and important in bone homeostasis [63, 64]. MKNK2 encodes the protein MNK2, which is a downstream kinase activated by MAPK1 [65]. All these genes are related to erythocytes but no common function emerges, which prohibits the identification of a core network gene, targeted by drugs to reduce simultaneously the level of SGE and the differentiation. Moreover, Artemisinin and Indomethacin decreased the noise of LDHA and MTFR1 genes, whereas MB-3 increased it. These genes are known to be involved in the metabolic switch that has been shown to be a key for the avian erythroid differentiation process [31]. Drugs could control SGE and differentiation though the modulation of the metabolic pathway needed to progress during the erythropoiesis. For further analyses, it could be interesting to further investigate these metabolic genes and the influence of their SGE change on the differentiation process.

It is important to note that our work focuses on genes that encode for erythroid differentiation and might not represent all genes. Overall it is advisable to use more caution when interpreting the importance of the role of gene affected by drugs in this study.

Thus, we can not exclude that there may exist another set of genes preferentially affected by drug that increase SGE.

We then investigated a potential role for variation in the mean gene expression that could explain the variation of the level of SGE.

Modifying the level of SGE is accompanied by a variation in the mean gene expression level for two drugs out of three. The decrease of mean gene expression under MB-3 treatment has been shown not to be significant in a different system [13]. Also, it has not been reported that Artemisinin affects mean gene expression in any other cellular system. However, the fact that Indomethacin treatment decreased gene-wise entropy clearly without affecting the mean gene-wise expression level reinforces the fact that the modification on the differentiation process is not associated to a modification in mean gene expression but only to a non-specific modulation of SGE.

Collectively, these results suggest that neither common genes nor common mechanisms could explain the observed effect of the three drugs simultaneously. This reinforces the fact that modulation of cell-to-cell variability is strongly accompanied by a change of differentiation, independently of a gene function or a specific mechanism involved.

This could be explained by adopting a dynamical systems view on the differentiation process, in the wake of Waddington’s proposal [66]. In such a view, we could consider that in the highly dimensional gene expression space, an equilibrium cell state could be compared to a valley in an epigenetic landscape [26]. It has been shown that entropy is a useful tool for analyzing stochastic processes [67] and distinguishes between equilibrium and transition states [50]. When we reduce SGE using Indomethacin or Artemisinin, we dig the valley, limiting the ability of cells to escape from a self-renewal equilibrium. Their probability to reach the new equilibrium state is reduced. Inversely, when we increase SGE using MB-3, we flatten the valley and improve the ability of cells to explore a larger dynamical landscape, and increase their probability to reach the new differentiated equilibrium state. Alternatively, we could think that drugs modulate differentiation dynamics, flattening or digging valleys, allowing cells to increase or decrease their probability to escape from the valley. Cells will experience a modulation of the amount of SGE as the consequence of their stability in the high dimensional gene expression space. Once cells achieve their journey, they stabilize their new gene expression pattern (the differentiated genetic profile) and return to a basal level of SGE [25, 26, 28]. In such a view, stochastic gene expression favours cells making the decision to differentiate, modifying the structure of the valley in which cells are moving. In a recent perspective, this same process of actively shaping the Waddington Landscape has been described in terms of a Plinko board, whose nail configuration, composition, and patterning can be modified towards forward stochastic design [12]. Similarly to our initial description [31], the variation of cell-to-cell gene expression in other differentiation systems has been recently described [13, 32–34, 68]. Furthermore, a strong evidence of the relation between transcriptional heterogeneity and cell fate transitions was demonstrated recently through destabilization of the histone acetylation leading simultaneously to an increase of SGE and differentiation of mouse embryonic stem cells [13] and myogenic progenitors [36]. This is fully backed up by our own data that also establish that the inverse (reducing simultaneously differentiation and SGE) can also be demonstrated.

Conclusion

We show in primary erythroid progenitor cells that a pharmacological modification of SGE is consistently accompanied by a modulation of the differentiation process. Similar experiments using the design principles described above can be used to help establish the contribution of variability to biological processes and to separate cause from consequence [45]. It could therefore be important to study the potential use of such drugs in differentiation-related diseases such as tumoral cell progression [69], as exemplified by chronic myeloid leukemia [12, 70, 71], paving the way to a “treatment by noise” of at least some cancer-related diseases.

Supporting information

S1 Fig [pdf]
interaction analysis between drugs’ targets.

S2 Fig [pdf]
Gene expression distributions.

S3 Fig [black]
The model reproduces the cellular kinetics observed .

S4 Fig [pdf]
Drugs do not change the erythroid commitment.

S5 Fig [pdf]
Model selection by Akaike’s weights.

S6 Fig [pdf]
Drug toxicity in T2EC.