A mathematical modelling framework for the regulation of intra-cellular OCT4 in human pluripotent stem cells
L E Wadkin, S Orozco-Fuentes, I Neganova, M Lako, N G Parker, A Shukurov
AA mathematical modelling framework for the regulation ofintra-cellular OCT4 in human pluripotent stem cells
L E Wadkin , S Orozco-Fuentes , I Neganova , M Lako , N G Parker , A Shukurov School of Mathematics, Statistics and Physics, Newcastle University, UK, NE1 7RU Department of Mathematics, Physics and Electrical Engineering, NorthumbriaUniversity, Newcastle upon Tyne, UK Institute of Cytology, RAS St Petersburg, Russia Bioscience Institute, Newcastle University, UK, NE1 3BZ* [email protected]
Abstract
Human pluripotent stem cells (hPSCs) have promising clinical applications inregenerative medicine, drug-discovery and personalised medicine due to their potentialto differentiate into all cell types, a property know as pluripotency. A deeperunderstanding of how pluripotency is regulated is required to assist in controllingpluripotency and differentiation trajectories experimentally. Mathematical modellingprovides a non-invasive tool through which to explore, characterise and replicate theregulation of pluripotency and the consequences on cell fate. Here we use experimentaldata of the expression of the pluripotency transcription factor OCT4 in a growing hPSCcolony to develop and evaluate mathematical models for temporal pluripotencyregulation. We consider fractional Brownian motion and the stochastic logistic equationand explore the effects of both additive and multiplicative noise. We illustrate the useof time-dependent carrying capacities and the introduction of Allee effects to thestochastic logistic equation to describe cell differentiation. This mathematicalframework for describing intra-cellular OCT4 regulation can be extended to othertranscription factors and developed into sophisticated predictive models.February 26, 2021 1/37 a r X i v : . [ q - b i o . S C ] F e b ntroduction Human pluripotent stem cells, hPSCs, have the ability to self-renew through repeateddivisions and to differentiate into a wide range of cell types, a property known aspluripotency. The pluripotency of hPSCs is their defining characteristic, central to theirtouted applications in drug discovery, regenerative and personalised medicine [1–6].However, hPSCs exhibit complex behaviour and the in-vitro control of theirdifferentiation trajectories is challenging.Pluripotency is controlled by an inter-regulatory network of pluripotencytransciption factors, PTFs, including the genes OCT4, SOX2 and NANOG [7–9]. Thedestabilisation of PTFs and their interaction with chemical signalling pathways result indifferentiation away from the pluripotent state and into a specialised cell [7, 10, 11]. Thisdecision of a cell to either remain pluripotent or to differentiate is known as its fatedecision. It is unknown how much cell fate decisions are led by inherited factors, asopposed to environmental factors and intra-cellular signalling as even clonal (geneticallyidentical) cells under apparently identical conditions make different fate decisions [12].In many in-vitro experiments the differentiation of hPSC populations is induced andfacilitated by a differentiation agent, such as BMP4 [13, 14].A narrow range of PTF expression is necessary to maintain cell pluripotency, withboth high and low expressions causing a shift from the pluripotent state [15, 16] andeven small fluctuations can bias cell fate decisions [17]. Furthermore, the PTFs areinherited asymmetrically as a cell divides, biasing the fate of the daughter cells andcontributing to colony heterogeneity [18–20] with the decision to differentiate largelydetermined before any differentiation stimulus is introduced [18]. Given the likely largenumber of factors involved in the fate decisions and our limited knowledge of theirnature, the probabilistic framework to modelling PTF dynamics appears to be the mostsuitable. However, careful, experiment-based quantification of the stochastic, temporaldynamics of PTFs is necessary to examine the resulting effects on cell fate.Statistical analysis and mathematical modelling are deepening our understanding ofhPSC behaviours and guiding the development of experimental protocols [21]. Recentmathematical models of cell pluripotency focus on describing the network of PTFs andthe resulting cell fate decisions to guide the optimisation and control of pluripotencyFebruary 26, 2021 2/37 n-vitro [21–23]. These models are informed by recent studies of fluctuations of PTFsthroughout colonies [17, 18, 24] and the spatial patterning of differentiation [25, 26].Many models use coupled differential equations based on the Hill equations [27]describing changes in concentrations of molecules to describe PTF fluctuations [28–30].Others use network analysis frameworks [31] or explore the mechanical aspects of thecell behaviour when both the model and data are complex [32]. These models often aimto describe the whole PTF regulatory network and it can be difficult to estimate themodel parameters accurately from experimental data [30].Here we focus on the methodology of building such mathematical models usingexperimental data for the transcription factor OCT4. Although the OCT4 dynamicswill be affected by many external factors and the remainder of the PTF network, thereare benefits to considering each PTF in isolation as the crucial first step; firstly, thissimplifies the model development process, allowing each element to be explored in asystematic way and secondly, the results provide a basis for comparison to the otherPTFs (e.g., NANOG and SOX2) from similar experiments. Similarly, althoughinteresting spatial patterning effects are seen in OCT4 [33], we will consider only theintra-cellular OCT4 behaviour through time. These simpler models can be used todescribe the stochastic nature of PTF regulation on shorter time scales and explore theeffects of each PTF on cell fate, before their development into coupled models of theentire pluripotency regulatory network.Here we systematically explore various mathematical models for the temporalregulation of the PTF OCT4. We aim to identify the optimal set of mathematical toolsrequired to reproduce the key quantitative features of experimental observations fromRef. [18] and the additional quantitative analysis of this dataset from Ref. [33]. Theframework discussed can be applied in future to other experimental datasets. Since PTFfluctuation is inherently stochastic [18, 24, 34, 35], we focus on different forms ofwell-established stochastic models to describe the behaviour, namely: fractionalBrownian motion and the stochastic logistic equation. The aim is to describe the PTFsas microstates before considering the macrostate of cellular pluripotency. Firstly, weintroduce the experimental data and outline the key features of OCT4 to be describedmathematically. Next, we explore fractional Brownian motion and the stochastic logisticequation for simulating temporal OCT4 before any cell differentiation occurs. WeFebruary 26, 2021 3/37onsider different types of random noise (additive and multiplicative [36, 37]) and theireffects. Finally, we examine the use of shifting carrying capacities and Allee effects tosimulate a reduction in OCT4 towards the differentiated state.
Experimental OCT4 fluctuations
We use experimental data of OCT4 expression in a growing hESC colony from Ref. [18]and our previous analysis of this data in Ref. [33] to guide model development.Although focused on one experiment, the mathematical framework outlined here iseasily adaptable to other experimental results. We use the experimental analysis inRef. [18] and Ref. [33] to illustrate the applicability of such models to PTF regulation.Here we summarise the experiment and main features of the data to be described by amathematical model.
Experiment summary
This experiment was carried out by Purvis Lab (University of North Carolina, School ofMedicine), and is published in Ref. [18]. The OCT4 levels (mean OCT4-mCherryfluorescence intensity) in a human embryonic stem cell colony were determined and cellswere live-imaged for 68 hours. The colony begins from 30 cells and grows over 68 hours(817 time frames) to 463 cells, with 1274 cell cycles elapsing within this time. After 40hours, the hESCs were treated with (100 ng/ml) bone-morphogenetic protein 4 (BMP4)to induce their differentiation towards distinct cell fates. The cell IDs, ancestries andpositions, ( x ( t ), y ( t )), were extracted along with their OCT4 immuno-fluorescenceintensity values (reported in arbitrary fluorescence units, a.f.u.). The measurements ofthe OCT4 signal at 5 minute intervals, results in a set of evenly sampled discreteobservations for each cell, OCT4( t ), OCT4( t ) , ..., OCT4( t n ), where t and t n denotethe times of cell birth and division, respectively. The values of t n range from 0.25–30hours across the population.To classify the cells as either self-renewing (pluripotent) or differentiated, theexpression levels of CDX2 were quantified at 68 hours. The final cells were fate classifiedusing CDX2 and OCT4 as either belonging exclusively to a pluripotent or differentiatedstate. A remaining group of cells were classified in an unknown category. Using theseFebruary 26, 2021 4/37ates, the cell population was traced back in time, spanning multiple cell divisions, witheach earlier cell labelled according to this pro-fate. In this paper we consider only thepluripotent and differentiated fate groups. Note that for times pre-BMP4 (before 40hours), the fate classification is a pro-fate based on the fate of the cells descendants. Temporal OCT4 features
The OCT4 expression of (pro-)pluripotent and (pro-)differentiated cells for the wholeexperimental time (68 hours) is shown in Fig 1(a). At 40 hours the differentiation agentBMP4 is added, after which there is a decline in OCT4 expression in the(pro-)differentiated cells. The (pro-)pluripotent cells retain their OCT4 expression levels.The distribution of all OCT4 expressions pre-differentiation is shown in Fig 1(b), withtemporal distributions in Fig 1(c) and (d) for pluripotent and differentiated pro-fatecells respectively. A detailed analysis of the experimental data is provided in Ref. [33].We identify several key features to capture in model development, summarised below.For simplicity, and due to the distinct behavioural differences identified pre- andpost-differentiation, we first consider modelling the temporal behaviour pre-BMP4before moving on to the effect of cell differentiation. • Pre-differentiation1. The time series exhibit stochastic noise with the Hurst exponent 0.38 in both(pro-)pluripotent and (pro-)differentiated cells, shown in Fig 1(a) andcalculated in Ref. [33]. A Hurst exponent < . ig 1. Experimental OCT4 properties. (a) The temporal OCT4 expression for all(pro-)pluripotent (purple) and (pro-)differentiated (green) cells up to 68 hours. At 40hours (dashed line) the differentiation agent BMP4 was added. Pre-differentiation: (b)The distribution of all OCT4 expressions for all (orange circles), pro-pluripotent (purplesquares) and pro-differentiated (green diamonds) cells. The distribution of OCT4expression for time intervals between zero and 40 hours for (c) pro-pluripotent and (d)pro-differentiated cells from the experiment. The colour bar shows the time of the bincentre. • Post-differentiation1. At the end of the experiment differentiated cells are classified according totheir OCT4 and CDX2 expressions. These differentiated cells shown apronounced reduction in OCT4 upon BMP4 addition (40 hours).2. There is a clear and natural separation between the two classified groupspost-BMP4 based on their OCT4 levels, with differentiated cells showingreduced OCT4 and pluripotent cells retaining OCT4 expression.
Results
Modelling OCT4 pre-differentiation
In the following sections we systematically explore the use of different stochastic modelsas a framework for temporal OCT4 regulation, aiming to capture the experimentalFebruary 26, 2021 6/37ehaviour described above and shown in Fig 1. All the models discussed have the samebasis, with the initial conditions and cellular division incorporated using the algorithmicbase model detailed below.
Base Model
1. We begin with a chosen initial number of cells, N = N , to match theexperimental conditions.2. Each of the N cells are allocated an initial OCT4 value. This is extractedprobabilistically from the kernel density fitting to the experimental distribution ofinitial OCT4, OCT4( t = 0), shown in Fig 2(a).3. Each of the N cells are allocated a cell cycle duration. This is extractedprobabilistically from the kernel density fitting to the experimental distribution ofcell cycle times for all pre-BMP4 cells, shown in Fig 2(b). Each cell’s startingposition in its cell cycle is chosen uniformly.4. For each of the N cells the OCT4 values for the duration of their cell cycle aresimulated using one of the stochastic models.5. Each of the N cells divide into two cells at the end of their cell cycle. For each ofthe two daughter cells, their initial OCT4 value is set to the pre-division OCT4value of the mother cell.6. Steps 4 and 5 are repeated for the number of required division events.When the cell cycle times are generated in step 3 it is necessary to specify how muchof the cell cycle has already elapsed. If all cells begin at the start of their cell cycle atthe start of the simulation then divisions will be synchronised, visible as ‘steps’ in thenumber of cells over time, N ( t ), as shown in Fig 2(c). Avoiding this synchronisation bystarting cells at different points in their cell cycle gives a more accurate representationof colony size, as shown in Fig 2(d).Although here we have used the analysis of the experimental data to inform theinitial conditions and the cell cycle simulation, this is flexible and can easily be adaptedto other experimental results. The OCT4 regulation itself is captured in step 4 and isFebruary 26, 2021 7/37 ig 2. The initial conditions and resulting population dynamics for thecommon base model. (a) The distribution of experimental initial OCT4 valueshistogram, OCT4( t = 0), with kernel density fitting shown in orange. (b) Theexperimental distribution of cell cycle duration times histogram for all cells pre-BMP4addition with kernel density fitting shown in orange. The number of cells over time, N ( t ), when cellular division is (c) synchronised and (d) not synchronised in step 3 ofthe common base model. Blue solid lines show the simulated population sizes withstandard deviation error range shown in light-blue (calculated over five realisations) andorange dashed lines show the experimental population.open to many mathematical modelling techniques. In the next section we use theexperimental results from Ref. [18] and [33] to systematically build a stochastic modelusing fractional Brownian motion and the stochastic logistic equation. Anti-persistent OCT4 fluctuations
One possibility for a simple model of OCT4 fluctuation is to assume that the expressionfluctuates symmetrically with no preferred trends or correlations. Mathematically thiswould be descried by a Wiener process, analogous to the physical phenomenon ofBrownian motion in one dimension and the starting point for many random walkmodels. However, the analysis of experimental OCT4 expression described above and inRef. [33] has shown that the OCT4 evolution is anti-persistent, with an average Hurstexponent of H = 0 .
38. This signifies that increases in OCT4 are more likely to befollowed by decreases, and vice versa. The Hurst exponent H (cid:54) = 0 . t , B H ( t ), with an initial value B H (0)and time increments B H ( t − s ) is defined by B H ( t ) = B H (0) + 1Γ( H + 0 . (cid:90) −∞ (cid:2) ( t − s ) H − . − ( − s ) H − . (cid:3) dB ( s )+ 1Γ( H + 0 . (cid:90) t ( t − s ) H − . dB ( s ) , (1)where H is the Hurst exponent and Γ is the gamma function [38]. There are severalways to simulate fBm, either exact or approximate [39–41]. Here we use the Matlabfunction ffgn [42] which uses the circulant embedding technique for H < . H > . wfbm (available in the Wavelettoolbox) which uses a wavelet based approximate simulation method [45].We can use fBm to simulate OCT4 over time (step 4 of the base model) with ascaling parameter σ which controls the level of noise, i.e., σB H . Example realisations ofthe fractional noise, corresponding fBm functions, and simulated OCT4 for varying H are shown in S1 Fig to illustrate the effect of the Hurst exponent. The parameter σ isestimated from the experimental data (for all pre-BMP4 cells) as the standard deviationof ∆OCT4 = OCT4( t ) − OCT4( t − σ ≈
90. Each time series for OCT4can then be generated as OCT4( t = 0) + σB H .For simplicity, we first consider both cell fates together with N = 16 cells, made upof 14 pro-pluripotent and two pro-differentiated cells to correspond to the experimentaldata [18]. For cells in the experimental colony H = 0 .
38 [33]. A comparable simulationusing fBm with 16 initial cells, H = 0 .
38, and σ = 90 is shown in S2 Fig. Note thatalthough we simulate from a limited number of starting cells, the number of OCT4values generated over 40 hours due to the 5 minute increments and cellular division isapproximately 30000. It is clear from S2 Fig that this level of anti-persistent regulationfrom the Hurst exponent is not sufficient to keep the OCT4 expression within the rangeseen in the experiment.February 26, 2021 9/37ne possible method of limiting the range of OCT4 is to impose boundaryconditions, such as absorbing or reflecting. For absorbing boundary conditions once theOCT4 level reaches the boundary the cell is theoretically removed in some way from theexperiment and its OCT4 time series does not continue. There is no indication orbiological evidence of particularly high or low OCT4 expressions resulting in cell deathexperimentally [18, 33]. However, high or low OCT4 expressions do accompany celldifferentiation [17], so the removal of cells via the boundary condition could correspondto the differentiation of cells if we consider pluripotent cells only. The OCT4 simulationfor fBm with absorbing boundary conditions is shown in S2 Fig.Reflecting boundary conditions imply that when the OCT4 expression reaches theboundary, it is reflected back in the opposite direction. Biologically this corresponds toan additional regulatory effect, either internal to the cell or external in experimentalconditions; if the OCT4 level in a cell becomes too low, there is systematic regulation toincrease it (and vice versa). The simulation using fBm with reflecting boundaryconditions is shown in S2 Fig. Reflecting boundary conditions produce a result moresimilar to the experiment than absorbing boundary conditions since cells are notartificially removed, but it still creates a sharper distribution boundary than seenexperimentally. Additionally, although the boundary conditions somewhat artificiallyforce the OCT4 into the desired range, the spread of the overall expressions is not wellcaptured.This illustrates that the anti-persistence from the Hurst exponent alone is notsufficient to capture the OCT4 regulation seen in the experiment, even with boundaryconditions. The imposition of any boundary conditions also requires furtherinvestigation to elucidate their nature, positioning and the biological implications.However, we can still incorporate fBm noise into other models to generate theanti-persistence. In the next section we consider describing temporal OCT4 with thestochastic logistic equation and explore the regulatory effects of a limiting carryingcapacity. The stochastic logistic equation
In this section we explore the application of the stochastic logistic equation (SLE) tosimulating temporal OCT4 regulation. The logistic equation is a widely used model ofFebruary 26, 2021 10/37opulation dynamics characterized by the growth rate of the population and its optimalsize called the carrying capacity. We adapt the logistic equation to the experimentaldata available, using the model for OCT4 variation, rather than the traditionalpopulation size. Since fBm alone does not fully capture the regulatory behaviour ofOCT4, some additional effects are clearly important. We consider the SLE withadditive noise, multiplicative noise, and the effect of a time-dependent carrying capacity.For simplicity, we again consider the two cell fates together initially.There are several ways stochasticity can be introduced into the logistic equation, e.g.,additive noise, multiplicative noise, a noisy growth rate parameter r or carryingcapacity K . Both additive and multiplicative noise can be used to regulate geneexpression [36]. The most straightforward of these is additive noise which can beintroduced by adding a noise term to the net rate of change in the PTF. The SLE withadditive random scatter to describe OCT4, O , over time, t , is then dOdt = rO (cid:18) − OK (cid:19) + σ A ξ, (2)where ξ is the stochastic noise (e.g., Wiener/Brownian noise, or fBM noise) and σ A is ascaling parameter controlling the magnitude of the scatter.Parameters that appear in Eq (2) are estimated from the experimental data(pre-BMP4). In keeping with the anti-persistence, the noise ξ corresponds to fBm noisewith the Hurst exponent H = 0 .
38 and the scaling parameter is again the standarddeviation of ∆OCT4, σ A = 90. We can also estimate the carrying capacity as themedian of all the experimental OCT4 values, K = 1290. The OCT4 dynamics simulatedusing Eq (2) with r = 0 .
02 is illustrated in Fig 3(a) and (b). Although the regulatoryeffect of the carrying capacity works well to capture the upper bound of OCT4expression, an additional boundary condition at small values of OCT4 is still required(if the stochasticity gives rise to
O < dO/dt < O → −∞ ). Adistinguishing feature not captured by the model is the positive skew in the distributionof all occurring OCT4 values, shown in Fig 1(b) and overlaid in Fig 3(b). The modelpromotes tighter regulation above the carrying capacity than below it, resulting in fewerOCT4 expressions above the carrying capacity than seen experimentally. This suggeststhat the stochasticity has some dependence on the current state of the system.February 26, 2021 11/37 ig 3. Comparison of experimental and simulated OCT4 using the SLEwith either additive or multiplicative noise. (a) Simulated OCT4 expression(blue) using the SLE with additive noise, Eq (2), with 16 initial cells, r = 0 . K = 1290, σ A = 90 and fBM noise with H = 0 .
38, with an absorbing boundarycondition at zero. The experimental OCT4 is shown in orange. (b) The correspondinghistogram of simulated OCT4 expression using Eq (2) with the experimental distributionand estimated carrying capacity ( K = 1290) in orange. (c) Simulated OCT4 expressionusing the SLE with multiplicative noise, Eq (3), with 16 initial cells, r = 0 . K = 1290, σ M = 0 . H = 0 .
38. (d) The correspondinghistogram of simulated OCT4 expression with the experimental distribution in orange.Whereas the additive noise in Eq (2) has no dependence on the state of the systemand corresponds to making dO/dt symmetrically noisy, multiplicative noise changesdepending on the current conditions. In the case of our temporal OCT4 simulation,multiplicative noise can be used to generate a scatter in the simulated data which has agreater magnitude when the system is close to the carrying capacity (thus resulting inmore stochastically high OCT4 expressions) and a reduced magnitude when far awayfrom the carrying capacity. Hints of this behaviour can be seen in Fig 1(a), with largerfluctuations apparent in the cells exhibiting above average OCT4 expression. Forsimulating the SLE with multiplicative noise we first consider the rearrangement of thelogistic equation, d ln( O ) dt = r (cid:18) − OK (cid:19) . February 26, 2021 12/37pplying the substitution X = ln( O ) and adding stochasticity ξ with noise scalingparameter σ M gives dXdt = r (cid:18) − e X K (cid:19) + σ M ξ, (3)which can then be used to simulate X = ln( O ), with the dynamics of OCT4 recoveredfrom O = e X . Example realisations of Eq (3) for both X and O are shown in S3 Fig toillustrate the effect of multiplicative noise in a typical logistic growth scenario forvarying σ M . The result is amplified noise for stochasticity occurring above the carryingcapacity.The temporal OCT4 dynamics simulated using the SLE with multiplicative noise,Eq (3), with fBM noise with H = 0 . r = 0 . K = 1290 and σ M = 0 . SLE with noise transition
Firstly, to capture the changing temporal skew for pluripotent cells, we could includeboth additive and multiplicative noise because different noise types reflect differentaspects in the cell behaviour [37] and both appear to be involved in the experimentallyobserved evolution of OCT4. If additive noise is dominant at early times, andmultiplicative noise at later times, the resulting OCT4 distribution will be symmetric atearly times and skewed at later times. We can consider the following rearrangement ofFebruary 26, 2021 13/37he stochastic logistic equation with additive noise d ln( O ) dt = r (cid:18) − OK (cid:19) + σ A O ξ , make the substitution X = ln( O ) and introduce the multiplicative noise term σ M ξ , dXdt = r (cid:18) − e X K (cid:19) + σ A e X ξ + σ M ξ . (4)As before, we can simulate the dynamics for X and recover the dynamics for O = e X .For simplicity, we can consider the change between additive and multiplicative noiseas a switch for pluripotent cells: for 0 < t <
20 h, σ A = 90 and σ M = 0, and for t >
20 h, σ A = 0 and σ M = 0 .
05. The additional parameters are specified in Table 1. Sincedifferentiated cells show reduced OCT4 expression throughout, they are given a lowercarrying capacity. The results for the OCT4 dynamics within this regime are shown inFig 4. The reduced carrying capacity for differentiated cells results in their lowerexpression throughout, shown in Fig 4(a). The overall OCT4 expression distributions inFig 4(b) are well described. The temporal distributions in Fig 4(c) illustrate the effectof the noise switch in the pluripotent cells, with the appearance of a positive skew atlater times, while the expression of differentiated cells in Fig 4(d) remains symmetricalat later times.
Table 1.
Simulation parameters for the OCT4 expression for pluripotent anddifferentiated cells using the SLE with both multiplicative and additive noise, Eq (4). At20 hours the noise switches from additive to multiplicative noise in the pluripotent cells.Parameter t <
20 h t ≥
20 hPluripotent N r K σ A
90 0 σ M H N r K σ A σ M H ig 4. The dynamics of OCT4 simulated using the SLE with a switchbetween additive and multiplicative noise. (a) The OCT4 dynamics between zeroand 40 hours for 14 pro-pluripotent (purple) and two pro-differentiated (green) initialcells following the SLE with both additive and multiplicative noise, Eq (4), with theparameters specified in Table 1. For pro-pluripotent cells the noise changes fromadditive to multiplicative at 20 hours. (b) The distribution of all simulated OCT4values for pro-pluripotent (purple) and pro-differentiated (green) cells with thecorresponding experimental distributions overlaid. The temporal distributions for (c)pro-pluripotent and (d) pro-differentiated cells split by time intervals.Although this model captures the overall distribution and provides the desiredtemporal change in skew (which could be further smoothed with a more sophisticatedtime-dependent noise function), it does not result in a shift in the mode expression asdrastic as the one apparent in Fig 1(c). For this we consider implementing atime-dependent carrying capacity in the next section. SLE with time-dependent carrying capacity
To reproduce the significant shift in the mode for the pluripotent cells we can employ atime-dependent carrying capacity. We use the stochastic logistic equation for all cells,with both multiplicative and additive noise, as in Eq (4), and a carrying capacity whichvaries with time, dXdt = r (cid:18) − e X K ( t ) (cid:19) + σ A e X ξ + σ M ξ . (5)We can estimate the carrying capacity as the median OCT4 between zero and 25February 26, 2021 15/37 able 2. Simulation parameters for generating OCT4 expression for pro-pluripotentand pro-differentiated cells using the SLE with additive and multiplicative noise, and atime-dependent carrying capacity, Eq (5).Parameter t <
25 h t ≥
25 hPluripotent N r K σ A σ M H N r K σ A σ M H K p ≈ K d ≈ K ≡ K p = K d ≈ ig 5. The dynamics of OCT4 simulated using the SLE with atime-dependent carrying capacity. (a) The OCT4 dynamics between zero and 40hours for 14 pro-pluripotent (purple) and two pro-differentiated (green) initial cellsfollowing the SLE with both additive and multiplicative noise and a time-dependentcarrying capacity, Eq (5), with the parameters specified in Table 2. For pro-pluripotentcells the carrying capacity reduces at 20 hours, whilst the carrying capacity forpro-differentiated cells is constant. (b) The distribution of all simulated OCT4 valuesfor pro-pluripotent (purple) and pro-differentiated (green) cells with the correspondingexperimental distributions overlaid. The temporal distributions for (c) pro-pluripotentand (d) pro-differentiated cells split by time intervals. Simulating cell differentiation
In the previous section we considered modelling temporal OCT4 regulation before anydifferentiation stimulus (BMP4) is added, corresponding to the time interval0 < t <
40 h in the experimental colony [18, 33]. The addition of BMP4 causes asignificant reduction in OCT4 expression in the differentiated cells, shown in Fig 1(a).The mean OCT4, shown in Fig 6(a) also shows the clear reduction in differentiated cells.The median and mode experimental OCT4 are shown in S4 Fig. We explore twomethods of modelling this reduction in OCT4 as differentiation is induced. Firstly, weapply the SLE with a time-dependent carrying capacity as discussed previously, andsecondly, we consider the use of the SLE with an Allee effect. Although not seen in thisexperiment, it should be noted that high OCT4 values can also correspond to celldifferentiation [17].February 26, 2021 17/37 ig 6. The experimental and simulated dynamics of OCT4 upondifferentiation at 40 hours.
The (a) experimental (i) OCT4 and (ii) mean OCT4with time. The (b) simulated (i) OCT4 and (ii) mean OCT4 with time withdifferentiation induced at 40 hours using a time-dependent carrying capacity, Eq (5),with the parameters specified in Table 3. The (c) simulated (i) OCT4 and (ii) meanOCT4 wth time with differentiation induced at 40 hours by introducing an Allee effectterm to the SLE, Eq (7), with r = 0 . K = 1290, σ A = 35, σ M = 0 .
035 and A = 1000. Differentiation with a time-dependent carrying capacity
We previously employed the SLE with a time-dependent carrying capacity, Eq (5), tosimulate a moderate reduction in the average OCT4 expression post-25 hours, as shownin Fig 5. We could extend this technique to simulate the more drastic reduction inOCT4 seen when the differentiation stimulus is added.As before, we can estimate the carrying capacities for t <
25 h as K p ≈ K d ≈ t >
25 h we canFebruary 26, 2021 18/37imulate the modest reduction in OCT4 expression for the pluripotent cells with areduction of the carrying capacity to K p ≈ K d ≈
300 in the time interval t >
40 h corresponds to cell differentiation.These shifting carrying capacities, along with the other model parameters are given inTable 3. The dynamics under this regime are shown in Fig 6(b) and S4 Fig. Thetime-dependent carrying capacity leads to the reduction of OCT4 in the differentiatedcell group, well capturing the dynamics of the experiment.
Table 3.
Simulation parameters for the OCT4 expression of pluripotent anddifferentiated cells using the SLE with additive and multiplicative noise, and atime-dependent carrying capacity, Eq (5), to capture induced differentiation.Parameter 0 ≤ t <
25 h 25 ≤ t <
40 h 40 ≤ t <
68 hPluripotent N r K σ A σ M H N r K σ A σ M H Differentiation with an Allee effect
Another possible method of modelling induced differentiation is the SLE with ademographic Allee effect. Allee effects are traditionally used for modelling populationnumbers, with the effect inhibiting population growth at low densities as observed inFebruary 26, 2021 19/37oth animal and cell populations [46–48]. The deterministic logistic equation for OCT4expression O with this effect incorporated has the form dOdt = rN (cid:18) − OK (cid:19) (cid:18) O − AK (cid:19) , (6)where A is critical point at which the Allee effect occurs. Note that there are othermethods of simulating Allee effects through e.g., difference equations [49, 50] andLotka-Voltera models [51, 52]. Here we use the logistic equation for consistency with ourprevious modelling results.The effect of the Allee term in Eq (6) on both dO/dt and the OCT4 expression O for an example system is illustrated in S5 Fig. For a weak Allee effect, A < O ( t = 0),the rate of change dO/dt remains positive for O < K but is significantly suppressed.For a stonger Allee effect,
A > O ( t = 0), dO/dt is negative for O < K and results in theOCT4 expression declining to zero. It is this declining effect we can employ to simulatethe reduction in OCT4 expression for the differentiated cells. The Allee effect can beintroduced at a certain time point resulting in either continued suppressed growth or adecline to zero. Examples of ‘switching on’ both weak and strong Allee effects duringlogistic growth are shown in S6 Fig.For simulating OCT4 expression through the differentiation process with the SLE,we can switch on the Allee effect term at the time the differentiation agent is added(40 h). If the OCT4 expression is below A , then the Allee effect will be strong and theOCT4 will decline to zero. The stochasticity in the system will mean that only some ofthe cells will meet this condition, with others having an OCT4 expression greater than A , and therefore continuing with (suppressed) logistic growth. The stochasticity willalso result in this effect taking place at all times past 40 h, so the differentiation processwill happen at different times for different cells. The SLE for X = ln( O ) with additivefBm noise ξ and multiplicative fBm noise ξ is dXdt = r (cid:18) − e X K (cid:19) (cid:18) e X − AK (cid:19) + σ A e X ξ + σ M ξ , (7)where A is the Allee effect critical point.The OCT4 dynamics for 16 cells simulated with the SLE, Eq (4), for t <
40 h andFebruary 26, 2021 20/37he SLE with an Allee effect, Eq (7), for t ≥
40 h with r = 0 . K = 1290, σ A = 35, σ M = 0 .
035 and A = 1000 are shown in Fig 6(c) and S7 Fig. Here the fates of each cellare identified at the end of the simulation, with the cells whose OCT4 has reduced as aresult of the Allee effect classed as differentiated, and the cells whose OCT4 hasremained constant as pluripotent. The model captures the reduction of OCT4 in thedifferentiated subset of cells whilst keeping a remaining pluripotent cell population.However, the OCT4 in the pro-differentiated group pre-Allee effect is no lower than forthe pro-pluripotent cell group, unlike in the experimental results. Furthermore, anadditional model would be required to introduce differentiated cells with high OCT4values. Discussion
We have explored different modelling techniques for describing temporal OCT4regulation, guided by previous analysis of experimental OCT4 expression in a growinghESC colony [18, 33], particularly fractional Brownian motion and the stochastic logisticequation. A differentiation agent, BMP4, was added to the cells at 40 hours and resultsin the reduction of OCT4 expression in the differentiated cells. Although not seen here,it is also possible for high OCT4 expression to accompany cell differentiation [17].Pre-BMP4 we identified some key features including an anti-persistent stochasticity, andfor pluripotent cells a temporal skew and shifting mode in the distribution of all OCT4expressions. All the models discussed follow a common base model which sets up theinitial conditions and describes cell proliferation. When adjusted to produceunsynchronised cell divisions, the base model describes well the population growth withtime, shown in Fig 2(d). We then focus on different mathematical methods ofgenerating the temporal OCT4 expressions for the cell population within this basemodel. The simulated populations consist of 16 cells (with 14 pro-pluripotent and twopro-differentiated) resulting in approximately 30000 simulated OCT4 expressions. Wehave taken a systematic approach, gradually building complexity to illustrate themethodology of developing statistical models for biological systems.Firstly, we consider modelling the OCT4 dynamics pre-BMP4, i.e., for t <
40 hours.The analysis in Ref. [33] revealed that OCT4 values fluctuate stochastically withFebruary 26, 2021 21/37nti-persistence and a Hurst exponent of 0.38, suggesting the use of fractional Brownianmotion (fBm) [38]. There is also further experimental evidence that gene expressionsand transcription factor dynamics display fractal characteristics [53]. The use of fBm isparticularly common in financial modelling [54–56], but it has also been used todescribe diffusion within crowded fluids (such as the cytoplasm of cells) [57] and thekinetics of transcription factors [58]. The stochasticity from fBm results in a widerrange of OCT4 values at later times than seen experimentally (an effect which isexacerbated with time).The range of OCT4 can be controlled artificially with boundary conditions (eitherabsorbing or reflecting), but the overall distribution of all OCT4 values is not wellcaptured, shown in S2 Fig. It is also unclear whether these boundary conditions arebiologically appropriate as OCT4 expression is regulated by a complex range of factorsacross the transcriptional, post-transcriptional and epigenetic regulationlevels [7, 11, 59, 60]. Interestingly, mechanical limits to transcription have been shown tonaturally generate bounds to transcriptional noise [61]. A boundary condition at zerocorresponds to the fact that OCT4 expression never becomes negative with the upperboundary representing a maximum possible value. Furthermore, what is the biologicalimplication of the removal of cells through through absorbing boundaries or the recoveryof expression through reflecting boundaries? One possibility for absorbing boundariesfor pro-pluripotent cells is to represent differentiation happening at both the upper andlower boundary [17]. Although fBm alone is not sufficient to capture the experimentalbehaviour, it does (by design) capture the anti-persistence ( H = 0 .
38) and so in all latermodel iterations we use fBm noise to generate the stochasticity.A somewhat less artificial method of keeping the OCT4 values within range is to usethe stochastic logistic equation (SLE), which has a regulating parameter of the carryingcapacity, K , which represents the maximum amount of OCT4 that can be expressedwithin each individual cell. Note that this could be due to limits on the expression ofOCT4 due to other members of the regulatory network which cause its down-regulation.In our model, the stochasticity allows for some fluctuations above K . Similarly to theboundary conditions this maximum value depends on the complex inter-regulatorynetwork of OCT4, however, we estimate the value of the carrying capacity from theexperimental results as the median of all OCT4 values (taking into account theFebruary 26, 2021 22/37tochasticity allowing for O > K ).There are many sources of noise within the system, with internal noise resultingfrom stochastic chemical reactions (represented by additive noise) and external noiseoriginating from fluctuations in other cellular components that indirectly cause variationin transcription factor dynamics [37]. We consider both additive and multiplicativenoise, shown in Fig 3. The introduction of multiplicative noise creates larger fluctuationsabove the carrying capacity, qualitatively similar to those seen in the experiment. Thisresults in a distribution of all OCT4 values well matched to the experiment, with theslight positive skew being captured. Both additive and multiplicative noise can be usedto regulate gene expression, with multiplicative noise allowing small deviations intranscription rates to lead to large fluctuations in protein productions [36].A property not captured by the SLE with either additive or multiplicative noise isthe time-dependency of this positive skew. It occurs only at later times, and only inpluripotent cells, shown in the time-discretised distributions of OCT4 in Fig 1(c). Thistemporal skew can be captured by the SLE with both additive and multiplicative noise,with the type of noise time-dependent; additive noise at early times producessymmetrical distributions of OCT4, with multiplicative noise at later times producingskewed distributions, shown in Fig 4. Here we changed the noise function stepwise, butthis could be further smoothed using a more sophisticated time-dependent noisefunction.Another interesting property of the experimental OCT4 is the decline in expressionfor pluripotent cells post-25 hours, shown in Fig 1(c). We consider capturing thisbehaviour using the SLE with a time-dependent carrying capacity. Since this parameteris likely to depend on a large number of biological factors, it is not unreasonable toexpect that it may change with environmental conditions and experimental time. Weconsider the pluripotent and differentiated cells separately, each with a differentcarrying capacity, corresponding to the suggestion that the decision to differentiate isdetermined pre-differentiation stimulus [18]. The carrying capacity for both cell groupsis reduced at 25 hours, resulting in a decline in OCT4 expression, particularly for thepluripotent cell group with originally higher expression. Although this technique welldescribes the experimental results (shown in Fig 5), it requires multiple parameterswhich need to be elucidated from further experimental data.February 26, 2021 23/37e then consider modelling the OCT4 regulation for all times, including the declinein expression due to the addition of the differentiation stimulus. We extend thetime-dependent carrying capacity approach, reducing the carrying capacity further forthe differentiated cell group at 40 hours. This well captures the decline in OCT4 upondifferentiation, along with the more subtle decline in pluripotent cells, shown inFig 6(b). Here we have used a stepwise change in the parameter K , but this is easilyadjustable to other experimental results and more sophisticated functions could be usedto capture other trends. Similarly, a population of high OCT4 differentiated cells couldbe introduced with a corresponding increase in their carrying capacity. Thepro-differentiated cells are identified from the outset and although this is notbiologically unreasonable, with evidence that cell fate is determined pre-differentiationagent [18], the model itself does not produce the two fate groups which limits its futurecapacity to develop into a predictive model.A method of inducing differentiation which naturally produces the two fate groups isthe SLE with an Allee effect. Allee effects are well used across mathematicalbiology [46–48], but we are not aware of their application to pluripotency transcriptionfactor expression. The Allee effect results in a decline to zero for cells whose OCT4expression fluctuates below the critical point A . The stochasticity in the system meansthat this condition is met for only some of the cells, causing the formation of adifferentiated cell group with reducing or zero OCT4 and a pluripotent cell group withstable OCT4 expression at the carrying capacity, shown in Fig 6(c). This model islimited to describing low OCT4 differentiated cells as seen in this experiment and highOCT4 differentiation would need to be incorporated through another technique. Thismodel could be combined with a time-dependent carrying capacity to capture thedecline in expression in pluripotent cells.The models discussed here are of a purely descriptive nature, but outline a possibleframework for modelling the regulation of OCT4. We have explored systematically awide range of effects that might be able to reproduce rather fine details in theexperimentally observed dynamics of the OCT4 expression and identified an adequateand optimal combination of such effects. However, the resulting model may not beunique and other approaches may be viable. To justify any model of this kind and todevelop it into a prognostic tool for in-silico experimentation, it should be assessed andFebruary 26, 2021 24/37ompared with targeted experiments. With this caveat, we believe that the modeldeveloped can be used as a provisional prognostic tool and basis for furthermathematical model development. A summary of the models discussed and theexperimental properties they capture corresponding to the key features identified isgiven in Table 4. Further time-lapse experiments monitoring PTFs are needed toconfirm which of these properties are inherent to OCT4 expression, and how they varydepending on experimental conditions, and to provide more extensive benchmarking forthe modelling approaches and assumptions. It will be informative to apply the samequantitative framework to the other predominant transcription factors, SOX2 andNANOG. Their individual regulatory dynamics could then be compared using the keydescriptive parameters, and any systematic differences identified. This information willhelp build the picture of the wider PTF system with the dynamics of the PTFsconsidered as part of an inter-linked network. In general, this highlights the need forfurther temporal experimental data on PTF regulation to build upon this mathematicalframework and develop more sophisticated predictive models. These models of themicrostate of PTF regulation will help inform longer time-scale models of thepluripotent macrostate. Table 4.
A summary of the key features identified experimentally and the models used to describe each behaviour.Key features ModelPre-differentiation 1. Stochastic noise with Hurst exponent of 0.38. fBm, Eq 12. Pro-differentiated cells show reduced OCT4 throughout. Incorporated through initial conditions.3. Positive skew of all pro-pluripotent OCT4 expressions. SLE (multiplicative noise), Eq 3 and 44. Reduction in pro-pluripotent OCT4 post 25 hours. SLE ( K ( t )), Eq 5Post-differentiation 1. Reduction in OCT4 expression for some cells. SLE ( K ( t ) or Allee effect), Eq 5 or Eq 72. Separation into pluripotent and differentiated groups. SLE (Allee effect), Eq 7 Acknowledgments
LEW would like to acknowledge support from the London Mathematical Society (EarlyCareer Fellowship). ML acknowledges BBSRC UK (BB/I020209/1). IN acknowledgesthe grant from the Russian Government 641 Program for the recruitment of the leadingscientists into 641 Russian Institution of Higher Education 14.w03.31.0029 and RFFIproject grant number 20-015-00060.February 26, 2021 25/37 eferences
1. Ebert AD, Svendsen CN. Human stem cells and drug screening: opportunitiesand challenges. Nat Rev Drug Discov. 2010;9(5):367–372.2. Zhu Z, Huangfu D. Human pluripotent stem cells: an emerging model indevelopmental biology. Development. 2013;140(4):705–717.3. Avior Y, Sagi I, Benvenisty N. Pluripotent stem cells in disease modelling anddrug discovery. Nat Rev Mol Cell Biol. 2016;17(3):170–182.4. Ilic D, Ogilvie C. Concise Review: Human Embryonic Stem Cells - What HaveWe Done? What Are We Doing? Where Are We Going? Stem Cells.2017;35(1):17–25. doi:10.1002/stem.2450.5. Shroff G, Titus JD, Shroff R. A review of the emerging potential therapy forneurological disorders: human embryonic stem cell therapy. Am J Stem Cells.2017;6(1):1.6. Trounson A, DeWitt ND. Pluripotent stem cells progressing to the clinic. NatRev Mol Cell Biol. 2016;17(3):194.7. Li M, Izpisua Belmonte JC. Deconstructing the pluripotency gene regulatorynetwork. Nat Cell Biol. 2018;20(4):382–392.8. Boyer LA, Lee TI, Cole MF, Johnstone SE, Levine SS, Zucker JP, et al. Coretranscriptional regulatory circuitry in human embryonic stem cells. Cell.2005;122(6):947—956. doi:10.1016/j.cell.2005.08.020.9. Chambers I, Tomlinson SR. The transcriptional foundation of pluripotency.Development. 2009;136(14):2311–2322. doi:10.1242/dev.024398.10. Kumar RM, Cahan P, Shalek AK, Satija R, DaleyKeyser A, Li H, et al.Deconstructing transcriptional heterogeneity in pluripotent stem cells. Nature.2014;516(7529):56–61.11. Wang Z, Oron E, Nelson B, Razis S, Ivanova N. Distinct lineage specificationroles for NANOG, OCT4, and SOX2 in human embryonic stem cells. Cell StemCell. 2012;10(4):440 – 454.February 26, 2021 26/372. Symmons O, Raj A. What’s luck got to do with it: single Cells, multiple fates,and biological non-determinism. Mol Cell. 2016;62(5):788 – 802.doi:https://doi.org/10.1016/j.molcel.2016.05.023.13. Kee K, Gonsalves JM, Clark AT, Pera RAR. Bone morphogenetic proteinsinduce germ cell differentiation from human embryonic stem cells. Stem CellsDev. 2006;15(6):831–837.14. Xu R, Chen X, Li DS, Li R, Addicks GC, Glennon C, et al. BMP4 initiateshuman embryonic stem cell differentiation to trophoblast. Nat Biotechnol.2002;20(12):1261–1264.15. Niwa H, Miyazaki J, Smith AG. Quantitative expression of Oct-3/4 definesdifferentiation, dedifferentiation or self-renewal of ES cells. Nat Genet.2000;24(4):372–376.16. Kopp JL, Ormsbee BD, Desler M, Rizzino A. Small increases in the level of Sox2trigger the differentiation of mouse embryonic stem cells. Stem cells.2008;26(4):903–911.17. Strebinger D, Deluz C, Friman ET, Govindan S, Alber AB, Suter DM.Endogenous fluctuations of OCT4 and SOX2 bias pluripotent cell fate decisions.Mol Syst Biol. 2019;15(9):e9002. doi:10.15252/msb.20199002.18. Wolff SC, Kedziora KM, Dumitru R, Dungee CD, Zikry TM, Beltran AS, et al.Inheritance of OCT4 predetermines fate choice in human embryonic stem cells.Mol Syst Biol. 2018;14(9):e8140. doi:10.15252/msb.20178140.19. Skamagki M, Wicher KB, J A, Ganguly S, Zernicka-Goetz M. AsymmetricLocalization of CDX2 mRNA during the First Cell-Fate Decision in Early MouseDevelopment. Cell Rep. 2013;3(2):442–457.20. Tee WW, Reinberg D. Chromatin features and the epigenetic regulation ofpluripotency states in ESCs. Development. 2014;141(12):2376–2390.doi:10.1242/dev.096982.February 26, 2021 27/371. Wadkin LE, Orozco-Fuentes S, Neganova I, Lako M, Shukurov A, Parker NG.The recent advances in the mathematical modelling of human pluripotent stemcells. SN Applied Sciences. 2020;2(2):276.22. Herberg M, Roeder I. Computational modelling of embryonic stem-cell fatecontrol. Development. 2015;142(13):2250–2260. doi:10.1242/dev.116343.23. Pir P, Le Nov`ere N. Mathematical models of pluripotent stem cells: at the dawnof predictive regenerative medicine. In: Systems Medicine. Springer; 2016. p.331–350.24. Torres-Padilla ME, Chambers I. Transcription factor heterogeneity in pluripotentstem cells: a stochastic advantage. Development. 2014;141(11):2173–2181.doi:10.1242/dev.102624.25. Rosowski KA, Mertz AF, Norcross S, Dufresne ER, Horsley V. Edges of humanembryonic stem cell colonies display distinct mechanical properties anddifferentiation potential. Sci Rep. 2015;5:14218.26. Warmflash A, Sorre B, Etoc F, Siggia ED, Brivanlou AH. A method torecapitulate early embryonic spatial patterning in human embryonic stem cells.NM. 2014;11(8):847–854.27. Hill AV. The combinations of haemoglobin with oxygen and with carbonmonoxide. I. Biochem J. 1913;7(5):471–480.28. Glauche I, Herberg M, Roeder I. Nanog Variability and Pluripotency Regulationof Embryonic Stem Cells - Insights from a Mathematical Model Analysis. PLoSOne. 2010;5(6):1–12. doi:10.1371/journal.pone.0011238.29. Chickarmane V, Troein C, Nuber UA, Sauro HM, Peterson C. TranscriptionalDynamics of the Embryonic Stem Cell Switch. PLoS Comput Biol.2006;2(9):1–13. doi:10.1371/journal.pcbi.0020123.30. Akberdin IR, Omelyanchuk NA, Fadeev SI, Leskova NE, Oschepkova EA,Kazantsev FV, et al. Pluripotency gene network dynamics: System views fromparametric analysis. PLoS One. 2018;13(3):1–24.February 26, 2021 28/371. Xu H, Ang YS, Sevilla A, Lemischka IR, Ma’ayan A. Construction andValidation of a Regulatory Network for Pluripotency and Self-Renewal of MouseEmbryonic Stem Cells. PLoS Comput Biol. 2014;10(8):1–14.doi:10.1371/journal.pcbi.1003777.32. Auddya D, Roth BJ. A mathematical description of a growing cell colony basedon the mechanical bidomain model. J Phys D Appl Phys. 2017;50(10):105401.33. Wadkin LE, Orozco-Fuentes S, Neganova I, Lako M, Barrio RA, Baggaley AW,et al. OCT4 expression in human embryonic stem cells: spatio-temporaldynamics and fate transitions. Phys Biol. 2020;doi:10.1088/1478-3975/abd22b.34. MacArthur BD, Lemischka IR. Statistical mechanics of pluripotency. Cell.2013;154:484–489.35. Holmes WR, Reyes de Mochel NS, Wang Q, Du H, Peng T, Chiang M, et al.Gene Expression Noise Enhances Robust Organization of the Early MammalianBlastocyst. PLoS Comput Biol. 2017;13(1):1–23.doi:10.1371/journal.pcbi.1005320.36. Hasty J, Pradines J, Dolnik M, Collins JJ. Noise-based switches and amplifiersfor gene expression. Proc Natl Acad Sci USA. 2000;97(5):2075–20 80.doi:10.1073/pnas.040411297.37. Liu XM, Xie HZ, Liu LG, Li ZB. Effect of multiplicative and additive noise ongenetic transcriptional regulatory mechanism. Physica. 2009;388(4):392–398.38. Mandelbrot BB, Van Ness JW. Fractional Brownian motions, fractional noisesand applications. SIAM Rev. 1968;10(4):422–437.39. Dieker AB, Mandjes M. On spectral simulation of fractional Brownian motion.Probab Eng Inform Sc. 2003;17(3):417–434.40. Dieker T. Simulation of fractional Brownian motion. Masters Thesis:Department of Mathematical Sciences, University of Twente. 2004;.41. Yin ZM. New methods for simulation of fractional Brownian motion. J ComputPhys. 1996;127(1):66–72.February 26, 2021 29/372. Stoev S. Simulation of Fractional Gaussian Noise *EXACT*; 2020. Availablefrom: .43. Dietrich CR, Newsam GN. Fast and Exact Simulation of Stationary GaussianProcesses through Circulant Embedding of the Covariance Matrix. SIAM J SciComput. 1997;18(4):1088–1107. doi:10.1137/S1064827592240555.44. Lowen SB. Efficent generation of fractional Brownian motion for simulation ofinfrared focal-plane array calibration drift. Methodol Comput Appl.1999;1(4):445–456.45. Abry P, Sellan F. The wavelet-based synthesis for fractional Brownian motionproposed by F. Sellan and Y. Meyer: Remarks and fast implementation; 1996.46. Drake JM, Kramer AM. Allee effects. Nat Edu Knowledge. 2001;3(10):2.47. Gascoigne JC, Lipcius RN. Allee effects driven by predation. J Appl Ecol.2004;41(5):801–810.48. Johnson KE, Howard G, Mo W, Strasser MK, Lima EABF, Huang S, et al.Cancer cell population growth kinetics at low densities deviate from theexponential growth model and suggest an Allee effect. PLoS Biol.2019;17(8):1–29.49. Elaydi SN, Sacker RJ. Population models with Allee effect: a new model. J BiolDynam. 2010;4(4):397–408.50. Wang M, Kot M, Neubert MG. Integrodifference equations, Allee effects, andinvasions. J Math Biol. 2002;44(2):150–168.51. Zhou SR, Liu YF, Wang G. The stability of predator–prey systems subject to theAllee effects. Theor Popul Biol. 2005;67(1):23–31.52. Lin Q. Allee effect increasing the final density of the species subject to the Alleeeffect in a Lotka–Volterra commensal symbiosis model. Adv Differ Equ.2018;196(1):1–9.February 26, 2021 30/373. Ghorbani M, Jonckheere EA, Bogdan P. Gene Expression Is Not Random:Scaling, Long-Range Cross-Dependence, and Fractal Characteristics of GeneRegulatory Networks. Front Physiol. 2018;9:1446.54. Cheridito P. Arbitrage in fractional Brownian motion models. Finance Stoch.2003;7(4):533–553.55. Xiao WL, Zhang WG, Zhang XL, Wang YL. Pricing currency options in afractional Brownian motion with jumps. Econ Model. 2010;27(5):935–942.56. Bender C, Sottinen T, Valkeila E. Fractional processes as models in stochasticfinance. In: Advanced mathematical methods for finance. Springer; 2011. p.75–103.57. Ernst D, Hellmann M, K¨ohler J, Weiss M. Fractional Brownian motion incrowded fluids. Soft Matter. 2012;8(18):4886–4889.58. Woringer M, Izeddin I, Favard C, Berry H. Anomalous Subdiffusion in LivingCells: Bridging the Gap Between Experiments and Realistic Models ThroughCollaborative Challenges. Front Phys. 2020;8:134. doi:10.3389/fphy.2020.00134.59. Shi G, Ying J. Role of OCT4 in maintaining and regaining stem cell pluripotency.Stem Cell Res Ther. 2010;1(5):39. doi:https://doi.org/10.1186/scrt39.60. Babaie Y, Herwig R, Greber B, Brink TC, Wruck W, Groth D, et al. Analysis ofOct4-dependent transcriptional networks regulating self-renewal and pluripotencyin human embryonic stem cells. Stem cells. 2007;25(2):500–510.61. Sevier SA, Kessler DA, Levine H. Mechanical bounds to transcriptional noise.Proc Natl Acad Sci USA. 2016;113(49):13983–13988.February 26, 2021 31/37 upporting information
S1 Fig. The effect of the Hurst exponent.
Realisations of simulated noise infractional Brownian motion with (a) H = 0 . H = 0 . H = 0 . B H (0) = 0. (g-i) Simulation of OCT4 for 40 hours, with ten initial cells, andtemporal OCT4 determined by simulated realisations of σB H with σ = 90 and (g) H = 0 .
1, (h) H = 0 . H = 0 . S2 Fig. The effect of boundary conditions on simulated OCT4.
SimulatedOCT4 expression (blue) using fBm with 16 initial cells, σ = 90 and H = 0 .
38 with (a)no, (c) absorbing, and (e) reflecting boundary conditions at zero and 2500. Theexperimental data is overlaid in orange. The corresponding histograms for simulatedOCT4 (blue) with (b) no, (d) absorbing, and (f) reflecting boundaries. The kerneldensity fitting to the experimental distribution is shown in orange.
S3 Fig. The SLE with multiplicative noise.
Realisations of the dynamics of (a) X = log( O ) and (b) O = e X from Eq (3) with r = 0 . K = 100, ξ = W t and σ M = 0(blue), 0.025 (orange) and 0.075 (yellow). S4 Fig. The average experimental and simulated OCT4 expressions usinga time-dependent carrying capacity.
The (a,c) experimental and (b,d) simulatedmedian and mode OCT4. The dynamics are simulated using the SLE with additive andmultiplicative noise, and a time-dependent carrying capacity, Eq (5), with parametersspecified in Table 3.
S5 Fig. The deterministic logistic equation with a demographic Alleeeffect.
The deterministic logistic equation with an initial condition of O = 10, r = 0 . K = 50 and an Allee effect, Eq (6), for (a) dO/dt and (b) O with A = 1(orange) and A = 50 (green). The deterministic logistic equation with no Allee effect isshown in blue. S6 Fig. Switching on a demographic Allee effect.
The deterministic logisticequation with an initial population size of N = 10, r = 0 . K = 50 (blue). TheFebruary 26, 2021 32/37llee effect term in Eq (6) is introduced at t = 25 h with (a) A = 20 (orange) and A = 25 (green) and (b) A = 40 (orange) and A = 50 (green). The deterministic logisticgrowth with no Allee effect is shown as blue dashed. S7 Fig. The average experimental and simulated OCT4 expressions usingan Allee effect.
The (a,c) experimental and (b,d) simulated median and mode OCT4.The dynamics are simulated using the SLE with an Allee effect at 40 hours, Eq (7),with r = 0 . K = 1290, σ A = 35, σ M = 0 .
035 and A = 1000.February 26, 2021 33/37 ig S1 February 26, 2021 34/37 ig S2Fig S3
February 26, 2021 35/37 ig S4Fig S5Fig S6
February 26, 2021 36/37 ig S7ig S7