Model definition

A host of experimental data has been collected which links PR-A, PR-B, and inflammatory drivers in the pregnancy uterus [7, 8, 20–22, 27, 32]. We translated the principles of the PR-A:PR-B hypothesis into equations which could be used to mathematically explore the dynamics and consistency of the biological hypothesis of how the progesterone receptors interact with inflammation during pregnancy. In essence, the PR-A:PR-B hypothesis describes a standard competitive interaction between the pro-pregnancy actions of PR-B and the pro-labor actions of PR-A where the activity of PR-A is related to the level of inflammation in the myometrium. As such, we chose to consider only PR-B and inflammation and incorporated the effects of PR-A into the inflammatory terms of the model (Fig. 1 a).

The growth of PR-B and inflammation is modeled by the terms \(\hat {b}\hat {B}(1-\frac {\hat {B}}{B_{c}})\) and \(\hat {i}\hat {I}(1-\frac {\hat {I}}{I_{c}})\) respectively. These terms imply that the levels of PR-B and inflammation increase in the presence of PR-B, \(\hat {B}\), and inflammation, \(\hat {I}\), at rates \(\hat {B}\) and \(\hat {I}\) respectively. The terms \((1-\frac {\hat {B}}{B_{c}})\) and \((1-\frac {\hat {I}}{I_{c}})\), impose a maximum or critical value on the level of PR-B and inflammation. At any given time, there is a finite level of PR-B induced activation of the transcriptional machinery. This critical level of PR-B is represented by the parameter B _c . Analogously, there is some saturable level of inflammatory drivers active in a myometrial cell, represented by the parameter I _c . The growth terms \(\hat {b}\hat {B}\) and \(\hat {i}\hat {I}\) in the equations for PR-B and inflammation will themselves increase in size as the amount of PR-B and inflammation increase, but in a way that is limited by the critical values for PR-B and inflammation. If \(\hat {B} =B_{c}\) or \(\hat {I} = I_{c}\), the limiting terms in parenthesis equal zero which causes the growth term to equal zero.

The depletion of PR-B and inflammation is modeled by the terms after the negative sign, namely \(k_{1}\hat {B}\hat {I}\) and \(k_{2}\hat {B}\hat {I}\). Qualitatively, this means that the rate of depletion of PR-B is the product of \(\hat {B}\), \(\hat {I}\), and a rate constant k ₁. The depletion of inflammation follows the same behavior with a different rate constant k ₂. The value of k ₁ accounts for the relative amount of PR-B to repress inflammation and k ₂ accounts for the relative impact of inflammation on PR-B activity. While we know that the phenomenon of PR-B repression of inflammation occurs, the exact mechanism for this repressive activity is not well understood. By allowing for k ₁ and k ₂ to take on different values relative to one another, we are able to explore multiple possible models for PR-B repression of inflammation.

Model nondimensionalization and simplification

where \(b = \frac {\hat {b}}{I_{c}k_{1}}, i = \frac {\hat {i}}{I_{c}k_{1}},\) and \(k = \frac {B_{c}k_{2}}{I_{c}k_{1}}.\)

Determining parameter values from transcriptome data

We obtained data from two published studies which examined transcriptional changes in the myometrium (obtained at the time of cesarean section delivery) of women who were not in labor (NIL: closed and rigid cervix and no indication of uterine contractions) and in labor (IL: cervix dilated >4 cm and rhythmic contractions). One study [28] in which transcriptome analysis was performed by microarray technology (henceforth referred to as the microarray dataset or microarray data) comprised 3 myometrial samples from NIL women at term (>37 weeks gestation), 3 samples from IL women at term, and 3 samples from IL women undergoing preterm (<37 weeks gestation) cesarean section delivery. Some confusion has emerged since publication about whether one of the term IL samples was truly in labor. We excluded this sample for our analysis of this dataset and combined the preterm and term IL samples into one IL group. The other study [29] used RNA sequencing (henceforth referred to as the RNAseq dataset or RNAseq data) of myometrium from 5 NIL women and 5 IL women at term. Both datasets collected independent, non-paired samples, comprising a total cohort of 18 unique samples (10 IL, 8 NIL).

We infer the activity of PR-B and pro-inflammatory drivers with two PR-B responsive genes to serve as surrogates for PR-B, FOXO1A and FKBP5 [20, 30], and three pro-inflammatory genes to serve as surrogates for inflammation, IL- 1β, IL-6 and IL-8 [31]. Although there may be more genes associated with PR-B and inflammatory response in the myometrium, a fully curated gene set does not exist for inferring the full set of transcriptome changes induced by PR-B and inflammation. Therefore, we focussed on a targeted set of well studied myometrium specific genes for inferring the activity of PR-B and inflammation in the model. The normalized values of the data for these genes N _j , are calculated using the equation, \(N_{j} = \frac {G_{j}- m}{M-m},\) where G _j is the value of the gene for patient j, M is the maximum expression value for that gene across patients in the dataset, and m is the minimum value of that gene across patients in the dataset. One consequence of this normalization procedure is that the transcriptome data has been nondimensionalized enabling integration of the dimensionless trascriptome data with the dimensionless constants. In this formulation the values of the surrogate genes parameterize the dimensionless model for each patient.

The normalization equation bounds the values for the PR-B and inflammation surrogates from 0 to 1 and makes the natural choice of values for critical levels of PR-B and inflammation B _c =I _c =1. We then assign the values of b and i using the normalized dimensionless values for the PR-B and inflammatory surrogate genes respectively and the value of k corresponds to the strength of PR-B’s anti-inflammatory actions. The values of b and i are determined by the normalized intensity of one of the PR-B and one of the inflammatory surrogate genes from the mRNA expression studies [28, 29] and there are six possible combinations of values for b and i depending upon which inflammatory and PR-B surrogate genes are chosen. Thus, for a particular patient, the values of b and i are determined by the normalized expression intensity from the RNA-seq or microarray study for one pair PR-B and inflammatory surrogate genes.

Calculating the probability of labor for each patient

Next, we quantify the behavior of k in order to apply our model to patient data. To do this, we have to derive the steady state solutions for our model. These solutions are the values of B and I which cause \(\frac {dB}{d\tau }=\frac {dI}{d\tau }=0\) and correspond to a state where the system undergoes no change. There are four steady states, also known as equilibrium points. These occur when the ordered pair for PR-B and inflammation, (B ^∗,I ^∗), is equal to, (0, 0), (1, 0), (0, 1), and \(\left (\frac {i(1-b)}{k-ib},\frac {b(k-i)}{k-ib}\right)=\left (\frac {i^{2}b(1-k)-ik(1-b)}{k(bi-k)}, \frac {kb(i-1)}{bi-k}\right)\). Of these, (0,0) is the trivial equilibrium point where neither PR-B nor inflammation is present, (1,0) is the quiescent equilibrium where PR-B is maximal and there is no inflammation, and (0,1) is the laboring equilibrium where there is no PR-B and inflammation is maximal. The quiescent equilibrium correspond to a PR-B dominant state and laboring equilibrium corresponds to an inflammatory dominant state. The fourth equilibrium point, the intermediate equilibrium, exists only for certain values of b and i between the quiescent and laboring equilibrium (For full derivation of equilibrium points see Appendix of derrivations).

Since the values of b and i are determined by PR-B and inflammatory surrogate genes scaled from 0 to 1, these terms are bounded to that interval. Furthermore, since B and I are bounded by B _c =1 and I _c =1, B ^∗ and I ^∗ are bounded to the square domain with vertices (0,0), (0,1), (1,1), and (1,0) and area 1. So, the intermediate equilibrium point, in both forms, should satisfy the constraints 0≤B ^∗≤1 and 0≤I ^∗≤1. By considering how this constraint impacts both forms of the intermediate equilibrium we can derive a set of constraints for the values of k. In order to allow for the full range of values of i, we find that i and k satisfy 0≤i<k≤1.

In the limit in the case where i=k, the intermediate equilibrium point equals the quiescent equilibrium point (1,0). If we visualize this state in phase space (Fig. 1 b) we see that all the vectors point away from the quiescent equilibrium point toward the laboring equilibrium. In phase space, these vectors define trajectories that indicate how the system would evolve in time given a certain starting point. The set of vectors pointing toward the laboring equilibrium point is known as the basin of attraction for the laboring equilibrium point.

We compute a probability of labor equal to the area of the laboring equilibrium point’s basin of attraction divided by the area of the domain, which in our case is 1. The area of the basin changes as k, i, and b change. This probabilistic interpretation of a phase space is reasonable under the assumptions that all possible pairs of values of B and I in the domain occur with equal likelihood. While it is clear that in a physiological context there are values of B and I that are more likely at different time points in pregnancy, since this information is not readily available our assumption of uniformity enables us to compute the probabilities in the most agnostic way possible. Given more information about the distribution of B and I over the course of pregnancy, it could be possible to get a more precise estimate of the probability of labor, but such an extension is not possible at this t.

Each patient in the microarray and RNA-seq datasets has an expression value for each of the surrogate genes, FOXO1A, FKBP5, IL- 1β, IL-6 and IL-8. In the absence of proteomic data precisely quantifying the protein level activity of these genes in vivo the mRNA expression levels can be combined with the framework of our mathematical model to approximate the functional activity of these genes at the time of labor, i.e. using FOXO1A or FKBP5 for PR-B and IL- 1β, IL-6, or IL-8 for inflammation. We calculated a probability of labor for each patient in each dataset using all six possible combinations of surrogate genes (FOXO1A, IL- 1β), (FKBP5, IL- 1β), (FOXO1A, IL-6), (FKBP5, IL-6), (FOXO1A, IL-8), and (FKBP5, IL-8) where each surrogate was used to set the parameters b and i in our model. We will hereafter refer to a pair of surrogate genes as a predictor. A probability of labor was computed for each patient which corresponds to the size of the basin of attraction for the laboring equilibrium point given a predictor pair of surrogate genes for b and i.

Classifier construction and assessment

After normalization within platform cohort, half of the IL and NIL samples from the total cohort of 18 were randomly selected as a training dataset to construct a classifier from five IL and four NIL samples. Probabilities of labor were calculated for each of the samples using a pair of predictor genes, one PR-B responsive and one inflammatory responsive, to value the parameters b and i in the model. Two nonparametric 95 % confidence intervals were computed for the the probabilities of labor for the NIL and IL samples. These intervals about the medians of the NIL and IL samples constituted the NIL and IL classifiers. If the intervals did not separate, then we discarded that classifier.

These 105,840 model classifiers (17,640 training set combinations × 6 combinations of model predictor genes) were compared to two types of null classifiers. We constructed single and two-gene null classifiers for each of the 17,640 combinations of samples using the normalized expression values of the raw datasets. The single gene null classifiers were built by constructing 95 % nonparametric confidence intervals for the IL and NIL training samples on the normalized expression values for the individual genes. The two-gene null classifiers were constructed by defining a two dimensional confidence region for the two predictor genes for the NIL and IL training samples. Precision, recall, F-score, and CSR metrics were calculated for all 264,600 null classifiers (17,640 training set combinations x 15 one or two-gene null classifiers) and compared to the model classifiers to assess i) the relative performance of the single and two-gene null classifiers and ii) the performance of our model’s two-gene classifier to the null two-gene classifiers. A full workflow of this approach can be found in Fig. 2.

Model results

The model has four equilibrium solutions, i.e. values for (B,I) pairs such that \(\frac {dB}{d\tau } = \frac {dI}{d\tau } = 0\). At these values of B and I, namely (0, 0), (1, 0), (0, 1), and \(\left (\frac {i(1-b)}{k-ib},\frac {b(k-i)}{k-ib}\right)\), the levels of PR-B and inflammation will remain constant. Each of these solutions corresponds to a physiological condition in the myometrium and the solutions can be stable, where the trajectories in the phase space in the neighborhood of the equilibrium point toward it, or unstable, where the trajectories in the neighborhood of the equilibrium point away from it.

The values of the parameters in the model influence the size of the basin of attraction for the labor equilibrium point and alter the stability of the equilibrium points. A tipping point exists in our model, where the quiescent equilibrium point transitions from stable to unstable and the phase space becomes completely biased toward the laboring equilibrium point. This change in stability is called a bifurcation and linear stability analysis allows us to compute the exact parameter values which will cause the stability to change. The end result of this analysis is a quantitative prediction for the tipping point of labor, the conditions under which the myometrial cell is permanently in the laboring phenotype leading to uterine emptying.

The point (0,0) corresponds to a myometrial cell that is not expressing any PR-B and has no inflammation. The solution (1,0) is the quiescent equilibrium corresponding to a physiological state where PR-B is at its maximal level with no inflammation. The solution (0,1) is the laboring equilibrium where inflammation is maximized and no PR-B is present. The last solution corresponds to an intermediate point between quiescent and laboring where the myometrial cell can pass into either phenotype.

where β=2(b B ^∗+i I ^∗)+k B ^∗+I ^∗−(i+b) and γ=i b(1−2I ^∗−2B ^∗+4B ^∗ I ^∗)+b B ^∗(−k+2k B ^∗)+i I ^∗(2I ^∗−1)+I ^∗ B ^∗(k−1). More details on the derivation of the eigenvalues can be found in the Appendix of derrivations.

The formula for the eigenvalues of the intermediate equilibrium indicates that three bifurcations are possible as i and b change. Firstly, if i=1 and b≤1 or if k≤i, then the quiescent equilibrium has one negative and one zero eigenvalue. In this case, the intermediate equilibrium point has moved through the phase space to collide with the quiescent equilibrium (Fig. 3). When these two equilibrium points combine, the quiescent equilibrium point changes from stable, where all temporal trajectories in the neighborhood of the equilibrium pointing toward it, to unstable with all the trajectories pointing away from the equilibrium point. This first bifurcation, the collision of the intermediate equilibrium with the quiescent equilibrium, corresponds to the physiological condition when the myometrium moves from quiescent to laboring. This condition corresponds to a probability of labor equal to one.

The converse of this is the second bifurcation that occurs when the intermediate equilibrium point collides with the laboring equilibrium point and the probability of labor is zero. This occurs when b=1 and i≤1, and results in a bifurcation where the laboring equilibrium point changes from stable to unstable. Physiologically this could correspond to a therapeutic intervention that preserves quiescence and prevents labor. Simulating therapeutic modulation of biomarkers that determine b and i could provide insight into what genes to modulate and what effect size is required to preserve quiescence.

When k=b i, the intermediate equilibrium becomes a singularity and is non-physiological. The third bifurcation occurs when b and i are both zero and the intermediate equilibrium collides with the trivial equilibrium point (Fig. 3). This third bifurcation may have physiological significance in the transition of the myometrial cell from non-pregnant to pregnant as inflammation and PR-B transition from inactive to active in the pregnancy uterus. However, since our model is based upon the activity of PR-B and inflammatory drivers during pregnancy, this bifurcation, though interesting, is beyond the scope of the present investigation.

The dynamical systems model was designed to explore the functional interaction between the anti-inflammatory actions of progesterone mediated by PR-B and the effect of inflammatory load on the contractile state of the human pregnancy uterus. In addition we sought to identify the conditions that induce a bifurcation in the model similar to the one that occurs when the uterus transitions to the laboring state. The dimensionless version of our model simplifies this task by enabling us to identify how changes in the three dimensionless parameters, b, i, and k, influence the trajectory of a hypothetical pregnancy phase space. The underlying rationale was that bifurcations correspond to physiologically important events in the timeline of pregnancy representing uterine quiescence and its transition to the laboring state. The interaction between the dimensionless model parameters i and k appears to be the most significant for initiating the labor bifurcation. The meaning of this interaction is that as long as the repressive capacity of PR-B, k, is greater than the activation rate of inflammation, i, quiescence will be maintained. This finding supports the PR-A:PR-B hypothesis since it recapitulates the important role of PR-B-mediated anti-inflammatory activity and shows how interference of this function, possibly by PR-A, destabilizes quiescence in favor of labor. Though we fixed k in the present analysis, we hypothesize that the level of k relative to i and b may be reflective of the trans-repressive activity of PR-A on PR-B. The phase space dynamics of modulating k seem to support this with higher values of k producing lower probabilities of labor and may be an interesting avenue for further investigation.

Predictive modeling of parturition datasets

We next applied the dimensionless version of our model with k=1 in Eq. 3 to predict the onset of labor in a cohort of 18 patient samples from myometrial transcriptome studies. Given a value of (b,i) for a particular patient, the model’s phase space reflects the state of the pregnancy myometrium for that particular patient. Interpreting the phase space as a probability means that we have a metric for predicting the likelihood of the patient going into labor. It is now possible to test the predictive power of particular sets of biomarkers for predicting the laboring phenotype by assigning the values of b and i based on the values of molecular markers of PR-B activity and inflammation. We can also assess the gain in robustness and predictive power of labor classifiers by comparing the performance of our model’s 105,840 possible classifiers gainst the 264,600 possible one and two-gene null classifiers built with the normalized gene expression data.

Identifying the specific inflammatory drivers that induce labor is an important step in identifying upstream and downstream therapeutic targets to delay the onset of premature labor. Interestingly, each of the inflammatory surrogate genes we used to construct our model’s classifiers has a subtly different biological function in the pregnancy uterus. IL-8 functions as a chemokine, drawing neutrophils and macrophages to tissues where it is expressed [32]. The lack of phenotypic predictability by IL-8, even when paired with a PR-B surrogate, may suggest that IL-8 does not play an important role in the inflammatory process of labor. In contrast IL- 1β and IL-6 [32] performed well as classifiers in our model when paired with PR-B surrogates, suggesting that the inflammatory processes associated with these genes are more important for labor onset than those of IL-8. In particular, IL-6 is a cytokine that plays an important role in both the canonical and non canonical JAK-STAT signaling pathway [33, 34], a pathway that may integrate effects of circulating and local myometrial cytokines. Recent work examining IL-6 as a blood-based biomarker for labor [35] provides further evidence that understanding and modeling this cytokine in the myometrium could be key to elucidating the driver pathways of labor.

The small sample sizes in the transcriptome datasets (N=18) cause us to exercise caution in the data analysis. The nonparametric methods of constructing confidence intervals and testing the separation of the NIL and IL groups are less powerful than comparable parametric testing, but were ultimately more appropriate due to their resistance to outliers, non normality of the data, and small sample size. Our approach of bootstrapping a distribution of classifiers and performance metrics allowed us to overcome the limitations of our cohort size and systematically test 370,440 classifiers and assess the performance of our model in a variety of contexts.

The present investigation of the predictive power of our model was limited by the availability of in vivo temporal data for the activity of PR-B and inflammation in the context of the pregnancy myometrium. If there were bloodborne biomarkers that could be used to infer the activity of these quantities, then it would be possible to test the application of our model in a clinical context modeling and predicting the pregnancy trajectories of actual women by measuring two biomarkers, one for PR-B and one for inflammation. This work is the first step in that direction by using transcriptome data from the myometrium to train a predictive model and demonstrate that our model provides the mechanistic hypothesis and framework to increase the predictive power of gene expression data.

Model nondimensionalization and simplification

where the values of b and i are determined by the normalized dimensionless values for the PR-B and inflammatory surrogate genes respectively and the value of k corresponds to the strength of PR-B’s anti-inflammatory actions.

Calculating the probability of labor for each patient

Since the values of b and i are determined by PR-B and inflammatory surrogate genes scaled from 0 to 1, these terms are bounded to that interval. Furthermore, since B and I are bounded by B _c =1 and I _c =1, B ^∗ and I ^∗ are bounded to the square domain with vertices (0,0), (0,1), (1,1), and (1,0) and area 1. So, the intermediate equilibrium point, in both forms, should satisfy the constraints 0≤B ^∗≤1 and 0≤I ^∗≤1. By considering how this constraint impacts both forms of the intermediate equilibrium we can derive a set of constraints for the values of k. In order to allow for the full range of values of i, we find that i and k satisfy 0≤i<k≤1.

Characterizing the stability of the quiescent and laboring equilibria

where (0,0) corresponds to a myometrial cell that is not expressing any PR-B and has no inflammation. The solution (1,0) corresponds to a myometrial cell that does not express PR-B and has no inflammation. The solution (1,0) is the quiescent equilibrium corresponding to a physiological state where PR-B is at its maximal level with no inflammation. The solution (0,1) is the laboring equilibrium where inflammation is maximized and no PR-B is present. The last solution corresponds to an intermediate point between quiescent and laboring where the myometrial cell can pass into either phenotype.

The sign of λ ₁ and λ ₂ determine the type and stability of each equilibrium point. The trivial equilibrium, (0,0), has two positive, real eigenvalues for all values of b and i indicating it is always unstable and can be classified as a source node. Physiologically, a myometrial cell in this state is not exposed to inflammatory stimuli and is not expressing PR-B. This state cannot endure long and like the equilibrium point is unstable, all trajectories are pointing away from the equilibrium. The quiescent and laboring equilibria have two real, negative eigenvalues each and are thus stable sink nodes. Both of these are stable so long as both b<1 and i<1. The intermediate equilibrium point allows us to identify the tipping point as i and b change. This equilibrium has two real eigenvalues, one positive and one negative, thus the intermediate equilibrium is a semi-stable, saddle node.

The formula for the eigenvalues of the intermediate equilibrium indicates that three bifurcations are possible as i and b change. Firstly, if i=1 and b≤1 or if k≤i, then the quiescent equilibrium has one negative and one zero eigenvalue. In this case, the intermediate equilibrium point has moved through the phase space to collide with the quiescent equilibrium. When these two equilibrium points combine, the quiescent equilibrium point changes stability from stable, where all temporal trajectories in the neighborhood of the equilibrium pointing toward it, to unstable with all the trajectories pointing away from the equilibrium point. Similarly, if b=1 and i≤1 the intermediate equilibrium collides with the laboring equilibrium resulting in a bifurcation where the laboring equilibrium point changes from stable to unstable. When k=b i the intermediate equilibrium becomes a singularity and is non-physiological. The third bifurcation occurs when b and i equal 0 and the intermediate equilibrium collides with the trivial equilibrium point.