Recapitulating human ovarian aging using random walks

Data sourcing and definitions

We performed a literature review of reports of numbers of primordial/nongrowing follicles, and established parameters for our modeling approach based upon (i) the distribution(s) of measured numbers of PFs present in the human ovary over time (Wallace & Kelsey, 2010) and (ii) the degree of variability between women in terms of their cessation of ovarian function (Weinstein et al., 2003; Nelson, 2009; Simpson, 2008; Persani, Rossetti & Cacciatore, 2010; McKinlay, Bifano & McKinlay, 1985; McKinlay, Brambilla & Posner, 1992). ANM data used in this study were generated from published plots (McKinlay, Bifano & McKinlay, 1985; McKinlay, Brambilla & Posner, 1992) using the data extraction tool WebPlotDigitizer (Rohatgi, 2021). Although loss of ovarian function prior to the age of 40 is defined clinically as reflective of a pathological state, termed primary ovarian insufficiency (POI; Nelson (2009)), we included the possibility that random action could result in measurable numbers of women that exhaust their PF reserve by that time.

Mathematical model of fluctuating ISR activity and PFGA

For a single PF in a given woman, we model fluctuating ISR activity by a RW. PFGA occurs when the ISR activity of this PF crosses a threshold. Crossing that threshold results in a PF’s subtraction from the ovarian reserve.

To describe our mathematical model precisely, let X(t) denote the ISR activity of a single PF at time t ≥ 0. Here, time t is the age of the woman so that time t = 0 corresponds to her birth.

To model ISR fluctuations as in Llerena Cari et al. (2021) (see also Fig. 1), suppose that the ISR activity in this PF either increases or decreases by an amount $\mathrm{\Delta }x$ > 0 over a time step $\mathrm{\Delta }t$ > 0. Mathematically, this is expressed as

(1) $$\begin{array}{c}X((n+1)\mathrm{\Delta }t)=\{\begin{array}{cc}X(n\mathrm{\Delta }t)+\mathrm{\Delta }x& \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{b}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}p,\\ X(n\mathrm{\Delta }t)-\mathrm{\Delta }x& \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{b}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}1-p,\end{array}\end{array}$$where p is the probability that ISR activity increases. The lefthand side of Eq. (1) is the ISR activity after n + 1 time steps, which is given by the ISR activity after n time steps (i.e., X(n $\mathrm{\Delta }t$)) plus or minus the amount $\mathrm{\Delta }x$. This type of model is called a random walk because the value of X “walks” by taking steps of size $\mathrm{\Delta }x$, and the direction of each step (either up or down) is random (Berg, 1993).

If the steps $\mathrm{\Delta }t$ and $\mathrm{\Delta }x$ are small, then the discrete random walk in Eq. (1) is equivalent to the continuous random walk whose dynamics are described by the stochastic differential equation (Gardiner, 2009),

(2) $$dX=-Vdt+\sqrt{2D}dW,$$where W(t) is a standard Brownian motion and

(3) $$D={\displaystyle \frac{{(\mathrm{\Delta }x)}^2}{2\mathrm{\Delta }t},V={\displaystyle \frac{2\mathrm{\Delta }x}{\mathrm{\Delta }t}\left({\displaystyle \frac{1}{2}-p}\right).}}$$

The equivalence of the discrete model Eq. (1) and the continuous model Eq. (2) for small steps $\mathrm{\Delta }t$ and $\mathrm{\Delta }x$ is shown in the Appendix (Supplementary Information).

In words, Eq. (2) means that over an infinitesimal time dt, the ISR activity X changes by an amount dX equal to a deterministic amount −V dt plus a stochastic fluctuation $\sqrt{2D}dW$, where dW is normally distributed with mean zero and variance dt. The parameter D is called the “diffusivity” as it describes the size of the stochastic fluctuations, and V is called the “drift” as it describes the deterministic (i.e., nonrandom) motion. Biologically, the drift V represents the efficiency of cellular repair in our model. Notice that V > 0 if p < 1/2, which means that X tends to decrease.

We suppose that a PF begins to grow when its ISR activity drops below some “growth” threshold. We also suppose that a PF dies before it begins to grow if its ISR activity rises above some “death” threshold. Without loss of generality, we take the growth threshold to be at X = 0 and the death threshold at X = L > 0. Hence, a PF leaves the reserve at the first time $\tau$ such that X(τ) ∉ (0, L). Since paths of X are random, this reserve exit time τ is random.

Let N be the number of PFs in a given woman’s reserve at birth, which we refer to as her starting supply. Let F(t) denotes the number of PFs in a given woman’s reserve at time t (and thus F(0) = N). We assume that the reserve exit times for each of the N PFs are independent and identically distributed, and thus the expected number of PFs in the reserve at time t is

(4) $$\mathbb{E}[F(t)]=N\mathbb{P}(\tau >t),$$where $\mathbb{E}$ denotes averaging over the reserve exit times.

Random walk model recapitulates PF decay

PF Data from −0.25 years to 0.1 years in Wallace & Kelsey (2010) were selected in order to establish the starting supply distribution (Fig. S1). The median starting supply was found to be

(5) $$N=3.23\times {10}^5.$$

For this median starting supply value Eq. (5), we sought parameter values for the random walk model to make the expected PF decay curve in Eq. (4) fit the PF counts reported by Wallace & Kelsey (2010). Since we can always rescale the ISR activity, without loss of generality we set the initial ISR activity to unity, X(0) = 1. This parameter search then yielded the following values for the diffusivity and drift,

(6) $$D=0.004{\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}}^{-1},V=0.051{\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}}^{-1},$$and any value of the death threshold L ≥ 2. By taking the limit of a large death threshold (i.e., $L\to \mathrm{\infty }$), we obtain the following formula for the expected follicle decay curve,

(7) $$\mathbb{E}[F(t)]={\displaystyle \frac{N}{2}[1+erf\left({\displaystyle \frac{1-Vt}{\sqrt{4Dt}}}\right)-e^{V/D}(1-erf\left({\displaystyle \frac{1+Vt}{\sqrt{4Dt}}}\right))],}$$where $erf(z):={\displaystyle \frac{2}{\sqrt{\pi }}{\int }_0^z{e}^{-y^2}dy}$ is the Gauss error function. In the Appendix (Supplementary Information), we derive Eq. (7) and show that the difference between Eq. (7) and the expected follicle decay curve for any L ≥ 2 is negligible for D and V in Eq. (6).

The formula Eq. (7) with the median starting supply value N in Eq. (5) and D and V in Eq. (6) yields the solid blue curve in Fig. 2A. In particular, by tuning only two parameters (D and V in Eq. (6)), this mathematical mathematical model yields a follicle decay curve in Eq. (7) that closely fits the human follicle decay data reported by Wallace & Kelsey (2010).

In fact, the model curve in Eq. (7) fits the PF data nearly as well as the non-mechanistic curve used in Wallace & Kelsey (2010). More precisely, Wallace & Kelsey (2010) posited a flexible functional form (the so-called asymmetric double-Gaussian cumulative (ADC) curve), and then chose five free parameters in this phenomenological function so that the resulting curve fit this same PF count data, plus some prenatal PF counts. Here, “fit” is defined in terms of the sum of squared errors between the curve and the logarithm of the PF counts. By this measure of fit, our model curve in Eq. (7) is only 3% worse than this prior curve.

Random walk model recapitulates human ANM distribution

Human menopause occurs when the number of PFs in the ovarian reserve drops below a threshold of approximately 1,000 (Faddy & Gosden, 1995). Hence, we can use our mathematical model of PF decay to study menopause timing.

The solid blue PF decay curve in Fig. 2A crosses the 1,000 PF threshold (horizontal black line) at age 51 years. Hence, this curve corresponds to an ANM equal to 51 years. However, there is considerable population variability in the starting supply N of PFs at birth. All else being equal, women with starting supplies higher (respectively, lower) than the median starting supply will tend to reach menopause later (respectively, earlier) than age 51.

To understand how starting supply population variability translates into ANM population variability, we first characterize the starting supply population distribution. As our starting supply data, we use the 30 PF counts in Wallace & Kelsey (2010) taken from women within a few months of birth.

We show in the Appendix (Supplementary Information) that these data are well-described by a log-normal distribution with parameters

(8) $$\mu =12.686,\sigma =\mathrm{0.497.}$$

That is, we model the distribution of the starting supply as

(9) $$N=\exp (\mu +\sigma Z),$$where Z is a standard normal random variable and μ and σ are in Eq. (8). The median of this starting supply distribution in Eq. (9) is exp(μ) = 3.23 × 10⁵ given in Eq. (5). As pointed out above, the solid blue PF decay curve in Fig. 2A that starts from this median starting supply yields an ANM equal to 51 years. Interestingly, the median ANM across a population is also approximately 51 years (Weinstein et al., 2003; McKinlay, Bifano & McKinlay, 1985; McKinlay, Brambilla & Posner, 1992). Hence, this solid blue curve in Fig. 2A can be understood as describing the “median woman.” The dashed blue curves in Fig. 2A show the PF decay curves for different quantiles of the starting supply distribution in Eq. (9) (namely, the 1%, 5%, 25%, 75%, 95%, 99% quantiles).

Model ANM output is compared with available information about the human ANM distribution as follows. First, the ANM distribution across a population of women who only vary in their starting supply N (as in Fig. 1A) is considered. In particular, the blue dashed curve in Fig. 2C was generated by simulating the PF decay dynamics of 10⁴ women, where each woman begins with a starting supply N sampled independently from the log-normal distribution in Eq. (9), and the ANM for each of the 10⁴ simulated women is the time when their PF reserve drops below 1,000. The PF decay dynamics for each woman follow the RW model described above, with diffusivity D and drift V in Eq. (6). The blue diamonds in Fig. 2C is the ANM distribution reported by Weinstein et al. (2003) when only variable starting supply was present in the RW model and drift was fixed, model output (blue dashed line) and the empirical (Weinstein et al., 2003). ANM distribution (blue diamonds) were in close agreement. We emphasize that this agreement follows merely from combining the empirical starting supply distribution in Eq. (9) with our RW model, where the two free parameters in the RW model (D and V in Eq. (6)) were chosen to fit PF decay data.

The agreement between the model and the ANM data in Fig. 2C is compelling, but there are some caveats. First, the ANM data from McKinlay, Brambilla & Posner (1992) in Fig. 2C (black squares) is evidently more variable than the ANM data from Weinstein et al. (2003), (blue diamonds). Second, the only source of population variability in the model for the blue curve in Fig. 2C is in starting supply. That is, each of the 10⁴ simulated women has identical diffusivity and drift parameters D and V in Eq. (6). However, it is not plausible that all intrinsic (e.g., genetic, epigenetic) and extrinsic (environmental) conditions that determine PFGA timing are entirely identical between women.

To address these two issues, we introduced an additional source of population heterogeneity into our model by allowing the drift parameter V to vary between simulated women in addition to the starting supply N. There are many possible choices that we could make for the population distribution of V, but for concreteness and we allow V to vary between women according to

(10) $$V=\overline{V}(1+cY),$$where $\overline{V}=0.051\mathrm{y}\mathrm{e}\mathrm{a}{\mathrm{r}}^{-1}$ as in Eq. (6), c = 0.03, and Y is a standard normal random variable independent of N in Eq. (9). In words, Eq. (10) simply means that the drift parameter V for each woman is normally distributed with mean (and median) $\overline{V}$ and 3% coefficient of variation.

The solid blue PF decay curve in Fig. 2B is identical to the “median woman” solid blue curve in Fig. 2A, since it is for the median starting supply in Eq. (5) and the median of the drift distribution in Eq. (10). The dashed blue curves in Fig. 2B are for the same quantiles of the starting supply distribution as in Fig. 2A, but the drift parameters for each of these six curves are sampled from the distribution in Eq. (10). Notice that these dashed curves sometimes cross, which means that a woman with a large starting supply may at some point have less PFs in her reserve than a woman with a smaller starting supply, due to differences in their drift parameters. Biologically, this represents individual differences in the rate of follicle loss over time as might be influenced by genetics or environmental exposures.

The black dotted curve in Fig. 2C shows the ANM distribution resulting from a variable starting supply and a variable drift. This ANM distribution was generated from the model analogous to the blue curve in Fig. 2C, except the drift parameter for each simulated woman varied according to Eq. (10). This plot shows that by introducing population heterogeneity into the drift parameter, the model yields an ANM distribution in line with the McKinlay, Brambilla & Posner (1992) data. Histograms depicting the ANM distributions produced by each of the RW conditions are shown in Fig. S2.

RW model flexibility

The analysis above shows that a simple RW model can recapitulate some prominent features of ovarian aging seen in nature. However, ovarian aging is a complex, multi-faceted process involving a variety of dynamic components. The purpose of this section is to show how our RW framework can be adapted to model different aspects of ovarian biology.

For example, the RW modeling framework can be used to investigate the effects of an acute drop in PF number. Such acute drops are common effects of certain cancer treatments (Fleischer, Vollenhoven & Weston, 2011; Das et al., 2012; Dunlop & Anderson, 2015; Chan et al., 2015). In Fig. 3A, RW traces for 50 simulated subjects are shown when Subject-Variable drift (identical to that shown in Fig. 2B) is applied, but PF starting supply is fixed at the reported median. The impact of an acute drop in PF number at approximately simulated age 12 is shown (asterisk, after De Vos, Devroey & Fauser (2010)). After PF number is adjusted in this way, RW traces are again shown for 50 simulated subjects, and an acceleration in the times that the ANM threshold is crossed is apparent when compared to uninterrupted decay. The RW modeling framework can also be used to investigate the purported acceleration of PF loss during human ovarian aging (Faddy et al., 1992) thought to arise due to declining levels of Antimullerian hormone (AMH) late in the third decade. While definitive experimental evidence that AMH slows PFGA and ovarian aging in other mammals is available (Durlinger et al., 1999; Pankhurst, 2017), it is less clear that AMH has this same role during ovarian aging in women in vivo (Leidy, Godfrey & Sutherland, 1998).

Acceleration in PF loss as might be seen in response to declining AMH levels can be incorporated into the model by allowing the drift parameter to increase over time. To illustrate, we let the mean drift parameter in Eq. (10) change at age 38 according to

(11) $$\begin{array}{c}\overline{V}=(\begin{array}{cc}0.024\mathrm{y}\mathrm{e}\mathrm{a}{\mathrm{r}}^{-1}& \mathrm{b}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{g}\mathrm{e}38,\\ 0.033\mathrm{y}\mathrm{e}\mathrm{a}{\mathrm{r}}^{-1}& \mathrm{a}\mathrm{f}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{e}38.\end{array}\end{array}$$

Given that the acceleration in PF loss around this time is more likely to be continuous, as the decline in AMH levels can be seen to be continuous in this window of time (Kelsey et al., 2011), we also tested more continuous drift modifications, but the difference between those and the simple stepwise drift modification was negligible (not shown). Figure 3B shows example PF decay trajectories for women with different starting supplies (equal to the quantiles in Fig. 2A) and whose drift parameters varied in time according to Eq. (11) and varied between women according to Eq. (10).

When Subject-Variable drift was applied and stepwise acceleration occurred (Fig. 3B), an ANM distribution was produced that again greatly resembles the actual ANM distribution (Fig. S2). In addition to a direct visualization of the impact of accelerating drift as might be expected due to declining AMH levels as women approach their forties, this approach demonstrates the flexibility of the RW model as it can be modified in order to begin to address these questions.

Software for further investigation

As illustrated above, a variety of biological factors can be investigated using the RW modeling framework. We therefore developed user-friendly R code for RW simulations with detailed documentation so that we ourselves and other interested parties can modify conditions, and, can add additional conditions that might influence PF loss (code available in a public repository, see the Methods section). This code allows one to alter a variety of conditions along possible degrees of freedom (e.g., variables like starting supply, acute follicle loss as shown in Fig. 3A, drift acceleration as in Fig. 3B, etc.) to investigate how such variables impact ovarian aging. Users can modify these conditions and add additional variables, perhaps as informed by their own experimental results.

Advantages

The RW approach is quite simple. Our model in its simplest form was able to fit both cross-sectional follicle decay data and an ANM distribution by choosing only two free parameters. The pattern of PFGA can be considered a by-product consequence of physiologically fluctuating ISR activity relative to a growth threshold, without the need for complex signaling or sensing that dictates whether an individual PF begins to grow or stays dormant. It may be that a more complex mechanism(s) ultimately dictates this decision, but the RW process appears capable of giving rise to the natural pattern of follicle loss.

Importantly, this interpretation in no way contradicts existing data on genetic and environmental factors that influence the overall rate of PFGA. Evaluation of genetic model systems and large-scale Genome Wide Association Studies (GWAS) have revealed genes and pathways that influence the rate of ovarian aging and the ANM (Stolk et al., 2012; Day et al., 2015). For example, Antimüllerian hormone (AMH; Durlinger et al., 1999; Pankhurst, 2017), and Phosphatase and tensin homolog deleted on chromosome 10 (PTEN; Reddy et al., 2008; Jagarlamudi et al., 2009; Cheng et al., 2015) negatively regulate PFGA as seen in accelerated loss of the PF reserve in mouse knockouts. Loss of Tumor necrosis factor alpha (Tnfα; Cui et al. (2011)), or, Tumor necrosis factor receptor 2 (Tnfr2; Greenfeld et al. (2007)) has the opposite effect, where PFGA is significantly slowed and the duration of ovarian function is extended. Downstream NFκB signaling also has been implicated in ovarian aging (Stolk et al., 2012), and disruption of NFκB inhibitory proteins IκBα and IκBβ also significantly slows the rate of PFGA (Wright et al., 2020). In addition to key protein signals, micro RNAs (miRNAs) and longer noncoding RNA species (ncRNAs) have been shown to regulate PFGA (Yan et al., 2019; Zheng et al., 2019; Maalouf, Liu & Pate, 2016; Wang et al., 2016; Yang et al., 2013). In the case of ncRNAs let-7/H19, this occurs due to their regulation of AMH levels (Qin et al., 2019). It has not been clear, however, how these and other signaling pathways combine in their action in order to determine known patterns of PF loss (Erickson, 2008; Kallen, Polotsky & Johnson, 2018; Zhang & Liu, 2015; Reddy, Zheng & Liu, 2010; Ernst et al., 2017). We now hypothesize that it is the combined action of all such identified factors that ultimately influence ISR activity, and that this integration of stress, damage, and signaling events gives rise to patterns of ovarian aging.

The recapitulation of ovarian aging patterns by RWs also illuminates several features of the biological process that are fascinating but have been difficult to explain. First, it is well-known that growing (preantral) follicles are present during prepubertal life in the human ovary, but insight as to what controls their numbers has been lacking. The RW approach provides a reasonable explanation for the onset and continuation of PFGA prior to puberty, and also the characteristic “plateau” in PF numbers during this time. This is because, if the RW process begins around the time of birth, there is a delay before the first follicles can engage in RWs for enough time such that they can cross the PFGA threshold in appreciable numbers. As presented here, appreciable numbers of PFs begin to be available in the window of time corresponding to menarche onset. Next, the supply of growing follicles between puberty and late in the third decade of life is consistent. The RW model gives rise to this stable supply of PFs in this window of (simulated) time, given the application of drift. It is also well-known that PFs and occasionally, growing follicles are present in the postmenopausal ovary. The RW model favors the continued presence and decay of PFs after menopause, and given a reasonably stable rate of PFGA, the rare commitment to growth would still occur in these later years. Mechanisms that control the maximum ANM can also be considered in new ways. Because some simulated subjects have nearly an order of magnitude greater starting supply of PFs, one or more mechanism must be in place that limits the ANM at approximately 62 years. Following along the results presented here, either a high PF starting supply must correspond to accelerated follicle loss compared to lower PF starting supplies, or the loss of the PFGA-slowing action of factors like AMH leads to acceleration in PFGA such that women almost never have more than 1,000 PFs remaining by age 62. These possibilities can be modeled as alterations in drift.

Our provision of the formal mathematical treatment of continuous and discrete-step models, and also user-modifiable code for discrete-step modeling makes it possible for interested parties to both reproduce our results, and also to test the effects of modifying the models upon patterns of PFGA and ovarian aging. Users can “tune” their modeling conditions as we did in cases of starting supply, drift, and diffusivity, with additional conditions that may better reflect actual signaling within the ovary, including mechanisms that have yet to be identified. Especially in the case of the discrete step RW R code, tuning can be performed by non-mathematics professionals. A final advantage of such tuning is evident in the way that inter-subject variability relates to patterns of ovarian aging and population-level ANM. While our results support starting supply as a plausible central factor in determining how the ANM differs between women (consistent with Depmann et al. (2015)), modeling output is improved when both starting supply and Subject-Variable drift are applied. The RW approach can be refined (and compared to alternative approaches) as biological mechanisms are revealed and higher quality female reproductive aging datasets become available.

Limitations

Despite what we consider compelling advantages, our approach has several limitations. First, we must acknowledge that the decades-long time scale of the process makes the normal behavior of individual PFs very difficult to track for direct follicle-by-follicle comparison to modeling results. In terms of the RW approach, while it is true that we only needed two free parameters to fit available follicle decay data and the known human ANM distribution, we acknowledge that there are many possible parameters (e.g., degrees of freedom) that could be tuned in order to produce RWs that faithfully reproduce these natural patterns. We also acknowledge that there may be differences between primordial follicles that influence their individual likelihood of growth activation in a more deterministic fashion. However, the RW model tuning that we performed was justifiable in terms of known biological measurements (starting supply) and the action of biological signaling factors (as in the modeling of accelerating PFGA due to declining AMH as the stepwise alteration of RW drift). It is true that the mathematical treatment of those biological features was done with the defined goal in mind of matching patterns of ovarian aging. Even so, it is satisfying that the action of a RW can so closely match PF loss and the ANM distribution when constructed in a logical fashion and subject to so few model variables.

We also acknowledge limitations with the best-available datasets used for RW modeling. As the authors that compiled the dataset(s) of histological PF numbers over time noted, how ovarian specimens were procured may influence the evaluated pattern of PF decline (Wallace & Kelsey, 2010). It is reasonable to expect that ovaries collected during autopsy are mostly representative of a random sampling from the general population. However, especially at later ages, ovary removal during elective surgery is less likely to be representative of the general population. Until non-invasive and accurate methods of PF number estimation become available, we and others will remain limited to cross-sectional data of PF numbers such as these. There are also caveats related to available datasets for the human ANM distribution. These are related to imprecision due to patient self-reporting their timing of menopause (defined as a year without menses; Weinstein et al., 2003; McKinlay, Bifano & McKinlay, 1985; McKinlay, Brambilla & Posner, 1992), and also to the fixed, artificial threshold of 1,000 PFs (Faddy & Gosden, 1995) used to define menopause onset in simulated women.

Next, despite expression of core ISR machinery in pregranulosa and granulosa cells, we have limited direct experimental evidence that ISR activity varies relative to a threshold in the pregranulosa cells of human PFs in vivo. We have generated immunofluorescence data that shows variable levels of the core ISR regulatory factor P-eIF2α in mouse PF pregranulosa cells and oocytes (Hagen-Lillevik, submitted, preprint available at https://www.researchsquare.com/article/rs-1682172/v1). Even so, despite our prior detection of P-eIF2α in human nongrowing follicles (Llerena Cari et al., 2021), a study comparing P-eIF2α levels in nongrowing follicles in multiple replicate human ovary specimens has not yet been performed. While it is highly likely that regional differences in stress and damage occur similarly between the mouse and human ovary, it may be that mechanisms that respond to these dynamic local conditions differ between the species. As mentioned, stochastic ISR checkpoint resolution is our favored model of the switch that activates PFGA, but the proposed RW could be greatly influenced by signals that function separately, or parallel, to the ISR.

Final considerations

RWs appear to provide a useful “systems biology” (Kitano, 2002) framework for the modeling of human ovarian aging. We consider the following final points the major implications of our approach. First and foremost, the recapitulation of ovarian aging by RWs suggests that random action may be an evolutionary strategy used to ensure that a minimum duration of ovarian function occurs in almost all women. It may be that more strict, deterministic PFGA control mechanisms would be more susceptible to dysregulation, and that randomness is protective against catastrophic PF loss. Next, our main model variables, PF starting supply and RW drift, are likely to be determined both by genetic/inherited factors in individuals and by environmental exposures. Genes identified in loss-of-function experimental studies and GWAS approaches may influence starting supply and/or drift. In the case of drift, these would likely be genes involved in the resolution of stress and cellular damage as at least in our model, these would modulate the probability of checkpoint resolution and cell cycle entry. Similarly, exposure to factors like chemotherapeutic agents could be acting at the level of drift, as the acute cellular response to damage would be upregulation of repair factors that again should impact the probability of checkpoint resolution, cell cycle entry, or death. Future experiments may be able to address these questions directly.

In the current model we treat the decision to undergo PFGA as if a single unit enters the cell cycle. As more information becomes available, we can more precisely model the multiple pregranulosa cells in human PFs, and how interactions between pregranulosa cells and with the oocyte influence RW behavior. We do not suggest that individual PFs behave in unpredictable ways. Instead, PFs likely respond to their local conditions in regulated, nonrandom ways, as dictated by known signaling pathways. PFs are exposed to conditions that change over time within the ovary, and simple, random processes are shown here to closely recapitulate the pattern of PFGA over time. We can therefore consider the mammalian ovary as a non-trivial, emergent, self-organizing system, with limited adaptive capacity. The ovary’s critical functions as the transmitting organ of the female germline as well as its support of female health and well-being prior to menopause are clearly non-trivial. That the simple system can deliver a consistent supply of maturing follicles that meets reproductive and endocrinological needs for decades is suggestive of emergent self-organization. Last, the ovary’s adaptive capacity is limited by the requirement that eggs produced can support reproduction by that subject. Adaptation can thus occur that alters ovarian function, but not in a way that compromises reproductive potential. Diverse future approaches including direct experimentation and further model refinement can be used to investigate how the pattern of ovarian aging in individual subjects and across populations of mammals, including women, is influenced by random behavior of ovarian follicles.