β2-Adrenergic Receptor Genotype and Other Variables that Contribute to Labor Pain and Progress

• β2-Adrenoceptor genotype affects the incidence of preterm labor, but its effect on labor progress and pain is not known.

• In 150 nulliparous parturients there was a strong association between β2-adrenoceptor genotype and progress of labor, with Asian women (who are likely to express one of the genotypes) to have slower labor progress.

• Black and/or heavy women were also noted to have slower latent labor.

ALTHOUGH term labor and delivery occurs tens of millions of times a year, it is poorly understood. Labor that begins prematurely—or proceeds either too slowly or too rapidly—poses serious risks for mothers and profound risks to fetuses. Labor progress is highly variable among women. Studies that have looked at demographic influences on labor progress have found a few significant differences that explain a small degree of this variability.1–4Slow labor progress is a common cause of cesarean delivery, the incidence of which has been increasing dramatically in recent years.5The high rate of cesarean section has caused a strain on healthcare resources and has likely resulted in increased incidence of abnormal placentation and postpartum hemorrhage.6Our previous work4with a biexponential mathematical model of labor demonstrated several factors associated with slower labor, including higher maternal weight, Asian ethnicity, and neuraxial analgesia.

Uterine contractility is modulated by several factors, including endogenous catecholamines that activate β2-adrenergic receptors (β2ARs). These receptors have two common polymorphisms in the coding region, at codons 16 (p.Arg16Gly) and 27 (p.Gln27Glu). Both have been shown to influence receptor function in a clinically significant manner.7We hypothesized that β2AR genotype may influence labor course.

Maternal labor pain is influenced by labor progress and varies greatly among women. We have previously shown that faster rates of cervical dilation predict faster pain development, and that Asian ethnicity predicts slower labor pain development.4The μ-opioid receptor (OPRM1) is a target for endogenous and exogenous opioid agonists. The gene has a common polymorphism at nucleotide 118 (codon 40, p.Asn40Asp), which is thought to impact baseline pain sensitivity8,9and has been shown to affect the efficacy of fentanyl spinal analgesia in human labor.10We hypothesize that this polymorphism may have an impact on labor pain.

We prospectively enrolled a nulliparous cohort from a single private practice. Subjects underwent quantitative behavioral testing for baseline pain sensitivity and were asked to provide demographic information as well as a DNA sample so that factors that affect labor progress and pain could be identified. This study was designed to use previously described mathematical models4,11to test the hypotheses that polymorphisms in β2AR would be associated with different rates of labor progress and that the p.Asn40Asp polymorphism in the opioid receptor would be predictive of labor pain. Furthermore, we hypothesized that we would be able to predict the magnitude of labor pain with quantitative behavioral pain testing during the third trimester of pregnancy.

This prospective cohort study was approved by the Institutional Review Board of Columbia University Medical Center (New York, NY). One hundred and fifty parturients were recruited from a single private obstetrical practice (S.D., N.J., P.R.) during a third trimester office visit. Eligible for inclusion were healthy nulliparous parturients aged 18–45 yr (gestational age 24–36 weeks) who planned vaginal delivery. Patients with intrauterine growth restriction, fetal malformations, preeclampsia, and planned cesarean sections were excluded. All subjects provided written informed consent. At the enrollment visit, patients completed a demographic questionnaire, underwent quantitative pain sensitivity testing, and provided a blood sample for genotyping.

Pressure sensitivity was evaluated using electronic von Frey sensory testing with an analgesiometer (IITC Life Science, Inc., Woodland Hills, CA). A pressure transducer was attached to a sterile, blunt 1-mm plastic tip. Graded pressure was applied via  the plastic tip to the bony aspect of the subject's nondominant shin. The subject was instructed to indicate when the stimulus first became painful. The pressure at which pain was first reported was recorded as the subject's pressure threshold.

Heat and cold sensitivity were tested using a TSA-II NeuroSensory Analyzer (Medoc, Ltd., Ramat Yishai, Israel). A 30-mm temperature probe was applied to the thenar eminence of the nondominant hand. A computer-driven program increased or decreased probe temperature 1° per second from 30–50°C (heat threshold) and 30–0°C (cold threshold). The subject was instructed to click a mouse button with the other hand when the temperature became painful; this process immediately terminated thermal stimulation. The temperature at which the test was terminated was recorded as the subject's heat or cold threshold. For each subject, each sensory test was completed in triplicate. Results were averaged for each individual.

A 10-ml blood sample was obtained from each subject and stored for DNA extraction and analysis. Genomic DNA was extracted with a Wizard® Genomic DNA Purification Kit (Promega Corporation, Madison, WI) according to manufacturer instructions. The genotype at two known polymorphic sites on the β2AR (rs1042713 and rs1042714, corresponding to codons 16 and 27) and one site on the OPRM1 gene (rs1799971, codon 40) was determined at the University of Geneva in Switzerland, as previously described.12,10 

Demographic variables, including maternal age, height, and self-reported ethnicity, were recorded at enrollment. Self-reported ethnicity was queried according to categories derived from the United States Census 2000. Patients were asked to choose only one category from the following choices: Asian, black, Hispanic, white, and other. Maternal weight, gestational age, infant weight, and treatment variables were recorded at delivery.

Labor was managed by three attending obstetricians from the practice group (S.D., N.J., P.R.) with support from obstetrics/gynecology residents. Generally, labor was managed with a standard clinical protocol; cervical examination was conducted every 4 h during latent labor and every 2 h during active labor unless otherwise clinically indicated. Timed cervical examinations were recorded. Labor induction, oxytocin augmentation, and time of membrane rupture were recorded. The time, type, and amount of pain treatments were recorded. Mode of delivery (e.g ., normal spontaneous vaginal delivery, operative vaginal delivery with vacuum or forceps, or cesarean section during the first stage or after full dilation) was recorded. Because full dilation is required in our model for temporal alignment of labor progress data, patients who underwent cesarean section during the first stage of labor were not included in our analysis.

On study enrollment, patients were instructed in the use of a numeric pain rating scale (NRS); scores ranged from 0 (no pain) to 10 (worst pain imaginable). They were told that they would be asked for pain scores at peak contraction. Pain scores were recorded by patients' nurses on arrival in the labor and delivery suite, at every change in condition, and before analgesia initiation. Pain scores were not recorded or analyzed for the purposes of this study after the patients received analgesia.

Because of our previous findings4regarding differences in labor pain and progress based on ethnicity, demographic and obstetric characteristics were compared among self-identified ethnic groups. Incidence data are reported as a percentage and compared using chi-square tests (Stata 10.1; StataCorp LP, College Station, TX). Distribution of the three genotypes was assessed according to ethnic group because of previously reported ethnic differences in β2AR13and OPRM114genotype incidence. To identify ethnic differences in genotypic distribution, each ethnic group's characteristics were compared with those of the total cohort using chi-square tests. Continuous data are reported as median (interquartile range). Differences in continuous variables are compared with a Mann-Whitney U test for two group comparisons or a Kruskal-Wallis test for three or more groups (GraphPad InStat 3.06, GraphPad Software, San Diego, CA). A P  value of less than 0.05 was considered significant in all analyses.

Labor progress was analyzed using the biexponential model for labor progress previously derived and validated by Debiec et al.  4using NONMEM (Nonlinear Mixed Effects Modeling; GloboMax LLC, Ellicott City, MD) with PLT tools (PLT Soft, San Francisco, CA). This structural model requires estimation of only three parameters: a rate constant for latent labor, a rate constant for active labor, and cervical dilation at which transition between the latent and active phases occurs.

Each genetic, demographic, and treatment variable was tested independently to determine whether it was predictive of a significant difference in active labor rate constant, latent labor rate constant, or cervical dilation at transition. The most parsimonious combination was determined for each effect. The factors tested for potential inclusion in the final model were β2AR genotype; maternal ethnicity, age, height, and weight; gestational age and infant weight; time of membrane rupture and mode of delivery; and oxytocin treatment and neuraxial analgesia. The factors that statistically improved the model individually were further considered for inclusion in the final model in rank order of statistical significance. After correction for the next strongest effect, factors that significantly improved the model were maintained in the final model.

Because genotype was significantly correlated with ethnicity (e.g. , none of the self-reported Asians had β2AR27-GG genotype), two models were constructed so that ethnicity and β2AR genotype could be considered separately after correction for other factors.

The bias in the model describes whether the model generally overpredicts or underpredicts actual measured values. Bias was calculated as the median prediction error (MPE). Therefore, MPE = median ([measured − predicted]/predicted). Inaccuracy was defined as the median absolute prediction error (MAPE). Therefore, MAPE = absolute value (median [(measured − predicted)/predicted]). MAPE is a measure of how close the typical observation is to the predicted value.

Mixed effects models are extremely sensitive and have the potential to detect real effects that are statistically significant, but perhaps too small to be clinically relevant. The effect size that is clinically relevant may differ among clinicians and clinical situations. To provide clinicians with a way to assess the absolute size of differences, taking statistical confidence into account, we constructed cumulative probability graphs. The details of their construction are described in the appendix. In brief, the value of the variable being compared (e.g ., latent labor rate constant, active labor rate constant, transition point, or the effect of a continuous variable like maternal weight on the previous variables) is displayed on the x-axis whereas the cumulative probability of the true value being at least that value is the denoted on the y-axis (from 0 to 1). The most likely value is located at 0.5. Thus, the probability of any clinically significant difference can be determined visually.

The effects of genetic, demographic, and treatment variables on labor pain were assessed using a sigmoid model with NONMEM and PLT tools, as described previously.4,15The sigmoid labor pain model requires four variables: (1) NRS^MIN, an NRS score that represents pain in early labor, (2) NRS^MAX, which represents pain in late labor, (3) γ, a slope function, which represents the rate of increase in pain during labor, and the closely related but unintuitive (4) CD⁵⁰, which is the cervical dilation at which pain has increased to 50% of maximal pain.

Similar to the modeling process for labor progress, each genetic, demographic, and treatment variable was first tested independently to determine whether it was predictive of a significant difference in NRS^MIN, NRS^MAX, γ, and CD⁵⁰. The best model was determined for each factor. In addition to OPRM1 genotype, the same demographic and obstetric factors were tested in the pain and progress models, except that neuraxial analgesia was not included because pain scores were not recorded after analgesia.

We considered variability between patients (interindividual variability) on NRS^MINand NRS^MAXas additive factors (i.e. , an individual's NRS^MINand NRS^MAXcould be larger or smaller than the nominal [population] value by an added or subtracted amount). In addition, we used an exponential model for the interindividual variability in CD⁵⁰. No interindividual variability was modeled for the slope function (γ) because it could not be estimated. The residual error that does not represent true differences between individuals, but rather random measurement error, was additive. Confidence in each covariate included in the labor pain model was assessed with cumulative probability graphs as described for the labor progress model.

The study was powered to detect a difference induced by a genetic covariate of 2.00 NRS units with a variance of 2, similar to the magnitude of effect we found in previous studies.15Although the prevalence of p.Arg16Arg homozygous genotypes is approximately 20%, it is approximately 11% for p.Glu27Glu. In addition, the minor allele (Glu) at codon 27 has a prevalence of approximately 35% whereas the prevalence of the minor allele at 16 (Arg) is approximately 45%.16For OPRM1, allele frequency for Asn (the A nucleotide) is approximately 30% in Caucasians and 60% in Asians.17,18Therefore, we recruited 150 patients who had approximately 100 vaginal deliveries with analyzable data. This sample size provides 80% power to detect a difference of 2.00 NRS units with a variance of 2 at a P  value of less than 0.05.

Two hundred and fifteen women were screened for enrollment between September 28, 2006, and November 6, 2008. One hundred and fifty women were enrolled, and 103 women had labor that progressed to full dilation. The analyzed dataset thus contains data from 103 women.

The screening and enrollment patient flow diagram is shown in figure 1. Of 103 patients with analyzable data, 14 were Asian, 11 were black, 10 were Hispanic, 67 were white, and one reported being in the “other” category. According to ethnicity, there were no differences in maternal age, weight, or height; gestational age at delivery; or newborn weight (table 1). Patients received 4–7 (median 5) cervical examinations during the first stage of labor. Twenty-six percent of patients had induced labor; 57% had spontaneous membrane rupture. Patients provided 1–3 (median 2) pain scores before analgesia. Ninety-six percent of patients had neuraxial analgesia initiated at some point during the first stage of labor (table 2).

All women were genotyped for two polymorphisms in the β2AR gene, p.Arg16Gly and p.Gln27Glu. Sixty-five of the 103 patients who reached full dilation were also genotyped for the p.Asp40Asn polymorphism in the OPRM1. Because genotyping of OPRM1 was added after study initiation, consent was not requested from all participants. Therefore, genotyping for the OPRM1 gene is incomplete. All resulting genotypes were within the expected Hardy-Weinberg equilibrium within the population as a whole and within each ethnic group. However, fewer Asians than expected expressed Glu at β2AR27 (P < 0.01; table 3). Asian ethnicity thus functions as a proxy for expression of Gln at the 27 position of β2AR. Because both Asian ethnicity and β2AR27 genotype were individually predictive of the transition point between active and latent labor, separate alternative final models were created to allow consideration of the effect of both factors individually because they would likely be collinear and inject instability to any model containing both variables. As has been previously reported,10,19Asians were also more likely than patients of other ethnicities to express the minor allele (G) at OPRM1 (P < 0.001).

The full dataset for all patients observed in the study (Supplemental Digital Content 1, https://links.lww.com/ALN/A685), NONMEM control files (Supplemental Digital Content 2, https://links.lww.com/ALN/A686) for the final labor progress and labor pain models, and diagnostic graphics of the model (Supplemental Digital Content 3, https://links.lww.com/ALN/A687) are provided online.

The biexponential model for labor progress was unbiased (MPE = 0.0 cm; MAPE = 0.94 cm; fig. 2A). When assessed independently, maternal weight, ethnicity, genotype at β2AR27 position, and the presence of neuraxial analgesia were significant covariates (table 4).

Because all of the Asian patients expressed at least one C at the p-β2AR27 position, genotype and ethnicity could not be considered together in one model. Accordingly, two final models were created: one optimized for ethnicity and the other for genotype. The results of the final model optimized for genetic information are shown in figure 2B. The genetically optimized model showed that patients with β2AR27 genotype CC (p.Gln27Gln) transition to active labor at 3.92 (2.95–4.50) cm whereas patients with other genotypes at this gene (CG and GG) transition to active labor earlier, at 2.73 (1.85–3.60) cm. The model also shows that heavier parturients had slower latent labor; their latent rate labor constants were lower by 0.0100 (0.0017–0.0200) h⁻¹per pound beyond median maternal weight. The active rate constant was very slightly increased by 0.007 (0.002–0.015) h⁻¹for each maternal pound. The epidural time scale factor is 0.29 (0.17–0.42), which means that labor progress is slower in the setting of earlier neuraxial analgesia. Specifically, the expected duration of labor in a parturient in which neuraxial analgesia is not present is approximately one third of the time observed in the presence of neuraxial analgesia.

The final model remained unbiased (MPE = 0.00 cm) and was made more accurate (MAPE = 0.74 cm) by considering the β2AR genotype, maternal weight, and neuraxial analgesia. At baseline, our biexponential model of labor progress was off by 0.95 cm. A residual error of 0.45 cm can be attributed to random variability such as measurement error, rather than to true differences between individuals (individual post hoc  Bayesian fit). As such, consideration of β2AR genotype, maternal weight, and neuraxial analgesia predicted 42% of the difference in labor progress between individuals, with 58% of the variability between individuals remaining to be explained by other, perhaps genetic or environmental, factors.

The final model optimized for ethnicity is shown in figure 2C. The ethnicity optimized model showed that black women have a very slow rate of latent labor 0.01 (−0.13 to 0.07) h⁻¹(table 4) compared with non-black women who have a latent labor rate constant of 0.08 (0.02–0.16) h⁻¹. Asian patients transition to active labor later than other patients, at 5.2 (3.7–7.0) cm compared with 3.30 (2.52–4.20) cm. The model optimized for ethnicity also finds that neuraxial analgesia is associated with prolonged labor progress by 0.32 (0.23–0.47) cm. The model is unbiased (MPE = 0.00 cm) and more accurate (MAPE = 0.72 cm) than the model without covariates.

Figure 3shows the cumulative probability distributions for each parameter tested in the genetically optimized model. This method was used to determine 95% CIs for each variable. Figure 4shows the probability distribution for a range of values for each variable tested in the ethnicity optimized model.

Figure 5Ashows simulations for labor progress in the genetically optimized model. The x-axis in figure 5Ashows the progress of labor from 1 cm cervical dilation (time runs forward in a normal way). The nominal patient is a patient with β2AR CG or GG genotype who is of median weight without neuraxial analgesia. Figure 5Bshows simulations for labor progress for the ethnicity optimized model.

Timed cervical dilations used in the pain model were derived from the labor progress model. Figure 6Ashows the basic sigmoid labor pain model without covariates. The model was unbiased (MPE = 0.00 NRS units), with a median inaccuracy (MAPE = 1.57 NRS units). Figure 6Bshows the fit of the final, covariate-adjusted model that relates labor pain to predicted cervical dilation. The final model was very slightly biased (MPE = 0.08 NRS units) and slightly more accurate (MAPE = 1.48) than the model without covariates. Individual parameter values for the pain model are found in table 5.

In the final model, patients who eventually required instrumental vaginal delivery had significantly higher pain scores in early labor, 5.43 (2.20–9.80) compared with initial pain scores of 0.79 (0.27–1.33; P < 0.001) among women who had a normal spontaneous vaginal delivery. Cold sensitivity, as determined with quantitative sensory testing, was a significant predictor for labor pain. Patients who were more sensitive to cold reported more labor pain. For every degree higher than the median cold threshold value, the predicted initial pain score (NRS^MIN) was greater by 0.14 (0.05–0.24) NRS units. Thus, the patient with the highest cold threshold (23°C) would be predicted to have an initial pain score that is 3.22 NRS units higher than the nominal patient. OPRM1 genotype was not a significant predictor in the final model. In the pain model, the error attributable to truly random factors (post hoc  Bayesian estimate of error) is 0.83 NRS units. The median absolute error in the initial model was 1.57, leaving 0.74 NRS units to be explained by true differences between individuals. Cold sensitivity and mode of delivery were able to explain only 12% of the variability.

Figure 7shows the probability of distribution of each variable contributing to the final pain model.

Figure 8is a simulation of maternal pain levels throughout the first stage of labor.

We have evaluated labor pain and progress in a prospectively enrolled cohort. Patients who express the CC allele at β2AR27 transition to active labor later than patients with other genotypes and thus have longer labors. Separately considered, Asian ethnicity predicted significantly slower labor progress. The much greater likelihood of expressing the CC allele at β2AR27 found in Asians may provide a genetic explanation for our previous and current finding of slower labor in Asian parturients.4Greater maternal weight and neuraxial analgesia were also predictive of slower labor, and is in accord with our and others previous findings.4Genetic data for a single polymorphism in β2AR allowed us to explain a significantly greater proportion of the residual variability compared with that which we were able to explain with demographic and treatment variables alone.4The labor pain model was minimally informed by cold sensitivity and mode of delivery; both variables were predictive of pain in early labor (NRS^MIN).

Because there was collinearity between ethnicity and genotype, we constructed two models to separately evaluate the two factors. It is well known that Gln is more commonly expressed in blacks,20and even more so in Asians than in the white population.21,22Indeed, in our simulations for labor progress, the effect of Asian ethnicity is very similar to the effect of expressing CC at the β2AR27 position (fig. 5, A and B). In fact, 11 of 14 Asian parturients had genotype CC at β2AR27 in our sample (none were GG). It may be that the association of Asian ehtnicity with slow labor progress is the result of this difference in β2AR, a hypothesis that cannot be adequately tested with such a high correlation of CC genotype and Asian ethnicity. Alternatively, β2AR could simply be a marker for Asian ethnicity and the causative difference in the rate of labor progress results from some other mechanism. We think this hypothesis is unlikely for two reasons: (1) non-Asians carrying the CC genotype experience a time course of labor progress that is similar to that of the Asian ethnic group, (2) β2AR activation (by endogenous catecholamines or the synthetic terbutaline) is well known to reduce uterine contractility, so a possible physiologic explanation is available. Because the β2AR27 polymorphism is thought to relate to receptor desensitization, it is mechanistically appealing to reason that the common allele is less desensitized and thus active and available to slow labor. The slow latent phase in black parturients may be related to β2AR genotype or to maternal weight. However, this issue is not clear and must be addressed in a larger, more diverse population. It should also be noted that the actual, in vivo  effect of these β2AR genotypes and haplotypes is still somewhat controversial.

Other genes, besides β2AR, might contribute to the variability observed in labor progress that remains unexplained. For example, the receptor for oxytocin has many polymorphisms that are thought to be important in autism spectrum disorders. However, the role of oxytocin-receptor polymorphisms in the normal progress of labor and pharmacologic administration of oxytocin for augmentation of labor has not yet been identified. We hope that genetic differences in β2AR, and potentially other receptors that modulate uterine contractility, will come together to describe a haplotype for labor progress. Using our models, we can and will evaluate the importance of other gene candidates in labor pain and progress in the future.

A recent analysis of the National Collaborative Perinatal Project database examined labor patterns in a large parturient population that delivered between 1959–1964.23Researchers23used interval-censored regression to estimate the distribution of times for labor progression. Comparing a simulation of their findings for nulliparas (Jun Zhang, M.D., Ph.D.) to ours, it is clear that the rates of latent labor are strikingly similar. However, patients in our cohort transitioned to active labor at a smaller cervical dilation than was observed in the historical cohort (fig. 9). This difference could be the result of the many demographic and treatment differences that have occurred during the past 50 yr. However, the difference is unlikely to be a result of the different modeling paradigms because these data from Zhang et al .23can be assessed easily with our structural model—as can our data with their high-order polynomial equations—with similar findings (data not shown).

We were able to explain a disappointingly small amount (12%) of the variability in labor pain that exists among women using the covariates we studied. Quantitative sensory testing for sensitivity to pressure, heat, and cold was performed in the third trimester of pregnancy. Although cold sensitivity was a significant predictor of early labor pain in our model, sensitivity to heat and pressure were not found to be important. It is unknown whether the transient receptor potential (TRP) channel subtypes, TRPM8 and TRPA1, which have a role in cold sensation,24are expressed in the uterus or play a role in labor pain. However, TRPV1 receptors are up-regulated in the human lower uterine segment and cervix at term.25These receptors may form heterodimers with TRPA1 subunits in sensory neurons.26Performing immunocytochemistry to detect cold receptors in human pregnant uterus will be helpful to elucidate this issue.

Genetic association studies that address pain focus on the OPRM1 gene as one of the highest ranked candidate polymorphisms for pain studies.9The μ-opioid receptor single nucleotide polymorphism OPRM1 118A>G (dbSNP rsl1799971) is associated with a 0.8-fold decrease in pressure pain intensity.27However, this study was likely underpowered to detect a difference by OPRM1 genotype. In this study, OPRM1 genotyping was added after study initiation. Because subject consent for OPRM1 genotyping was not requested in all written informed consent documents, there is incomplete genotyping in our cohort. Only 65 subjects were genotyped for OPRM1 gene. If the phenotype is conferred by the minor allele, or by the homozygous minor allele genotype, a much larger dataset will be required to detect a significant effect.

Aside from OPRM1, there are other genes known to be important in baseline pain sensitivity that have been suggested as candidate genes for pain trials.8These include transient receptor potential cation channel, subfamily V, member 1 (TRPV1); transient receptor potential cation channel, subfamily A, member 1 (TRPA1); catechol-O-methyl transferase; fatty acid amide hydrolase; and δ-opioid receptor OPRD1. Among these, polymorphisms in TRP genes may be of particular interest to researchers because they have been associated with differences in cold sensitivity in previous human trials and are highly expressed in the uterine cervix at term pregnancy.

It should be noted that there were significant limitations to our study. Our database was relatively small—consisting of 103 patients who reached full dilation—and we only tested for few candidate genes. Our trial was limited to nulliparous women from one private practice. This study design allowed us to limit variability in obstetrics practice, but it also limited demographic variability. Ideally, future studies will be conducted in larger cohorts or DNA-linked databases to allow more thorough evaluation of genetic effects on labor pain and progress. Also, as discussed extensively in Debiec et al .,4our pain model suffers from both left and right censoring because women often come to the labor room because of pain and are more likely to request neuraxial analgesia earlier if they have more pain. The use of this model in a population in which neuraxial anesthesia is less prevalent will help assess factors that affect pain in later labor. It is noteworthy, however, that in two separate cohorts of laboring women, one retrospectively and one prospectively enrolled, we have identified Asian ethnicity (now with a potential mechanistic explanation), maternal weight, and neuraxial analgesia as being predictive of labor progress.

A benefit of our structural model is that, in the future, it could be used with concurrent patient data to predict the time of full dilation. Prospective validation of this technique will allow us potentially to differentiate the criteria that are associated with successful vaginal delivery from those associated with cesarean section for dystocia. A better understanding of these criteria, combined with a reasonable estimation of the time required to achieve full cervical dilation, may help practitioners and patients have a more informed discussion about delivery plans.

The goal of modeling is to identify the most likely model that describes a particular dataset. NONMEM does this by calculating how likely each data point is. The calculated value is called the likelihood . The likelihood is virtually identical to the probability of the observation, but with a subtle nuance. Probability  is the term used to describe how likely a particular observation is within a particular model. Likelihood  is the term used to describe how likely a particular model is, given an observation. Because the calculation is the same, the terms will be used interchangeably.

The probability of several independent events is the product of the probabilities of each event. For example, if you flip a coin, the chance of it being “heads” is 0.5. If you flip it twice, the chance of it being heads both times is 0.5 × 0.5. The same is true for likelihood. The likelihood for the model, given all of the data, is of the product of the likelihood of each observation. However, the likelihood of each observation is a number between 0 and 1. If there are 1,000 numbers, and the likelihood of each observation is less than 1 (by definition), then the final product will be much less than 1. For example, let's say that the likelihood of each observation is 0.5. In that case, if there are 1,000 observations, then the likelihood of the model is 0.5¹⁰⁰⁰or 9.3 × 10⁻³⁰²(i.e ., a very small number). Calculation errors are possible, even with modern computers handling these small numbers. Fortunately, there is an easier approach. If one takes the log of each individual likelihood, the numbers can be added instead of multiplied. This process allows the computer to deal with more reasonably sized (rather than incredibly puny) numbers. For example, the log of 0.5 is −0.693. Adding up 1,000 of these gives us −693, a manageable number. It is much easier for computers to deal with −693 than 9.3 × 10⁻³⁰², even though log (9.3 × 10⁻³⁰²) =−693.

The most likely model has the maximum likelihood, and, of course, the highest log likelihood. However, search algorithms typically seek the minimum value of the function, not the maximum value. Thus, programs that search for the most likely model typically do so by minimizing the negative log likelihood.

Instead of minimizing the negative likelihood, NONMEM minimizes −2 log likelihood. The reason is the likelihood ratio test, a common test for statistical significance when looking at models. Consider two models, one of which has an additional parameter. For example, one model does not include age but the other model has an additional parameter that allows incorporation of age. We can test the statistical significance of the additional parameter by looking at the ratio of their likelihoods, L¹/L², where L¹and L²are the log likelihood values for each model. It turns out that the distribution of −2 log (L¹/L²) follows a chi-square test distribution. We can use the chi-square test distribution to see if the parameter is statistically significant. We can expand −2 log (L¹/L²) to −2 (log [L¹]− log [L²]) or (−2 [log (L¹)]) − (−2 [log (L²)]). Because NONMEM reports out −2 log likelihood (e.g. , −2 [log (L)]), −2 log (L¹/L²) is simply the difference between the NONMEM objective functions for the two models. This difference follows a chi-square test distribution, with the degrees of freedom being the difference in the number of parameters. For example, with a probability of 0.05, and one degree of freedom, the value of the chi-square test distribution is 3.84. Thus, if the difference in −2LL values (i.e. , the difference in NONMEM objective function values) for two models that differ by only one parameter exceeds 3.84, the parameter is significant with a P  value of less than 0.05.

The likelihood profile is a graph of the −2 log likelihood of the model, as a function of the value of the parameter. Typically, it is parabolic in shape, with the nadir of the parabola representing, by definition, the most likely value of the parameter. Figure 10shows a representative likelihood profile. The x-axis is the value of the parameter; the y-axis is −2 log likelihood. The best model, meaning the most likely model, is the model where the parameter equals 5, because that is the nadir of the parabola.

We can use the likelihood profile to determine how likely any possible value of the model is. One way is to simply subtract the nadir value of the log likelihood from the curve, so that it just touches the Y = 0 line at the nadir. That curve is shown in figure 11. In figure 11, the likelihood parabola intersects with the x-axis at the optimal value, in this case 5. The y-axis shows the number of log likelihood units that every point varies from the optimal value. A blue horizontal line is drawn at 3.84. The value of the chi-square test distribution for is a P  value of 0.05 with one degree of freedom. The thin black vertical lines are at the intersection of the log likelihood profile with the value 3.84, representing parameter values of −3.77 and 13.77, respectively. We can translate this result directly into a CI; there is a 95% probability that the value of the parameter is between −3.77 and 13.77. The remaining 5% is divided between values above and values below the 95% CI, as shown.

In a typical analysis, the question considered is whether the parameter value is statistically different from 0. We can see that 0 is within the 95% CI in figure 11, so the best value of the parameter, 5, is not statistically different from 0. However, let's assume that the parameter is of grave clinical importance, such as the number of years of additional life in a model that compares two chemotherapy regimens. Knowing the 95% CI, or the fact that the CI includes 0, does not provide the information we are looking for. We want to look at different possibilities for survival differences, how likely they are, and weigh that information when deciding which chemotherapy regimen to recommend.

To do this, we can convert the log likelihood profile in figure 11into a cumulative probability distribution. For values below the nadir, this number is half of the P  value associated with the −2 log likelihood difference. For values above the nadir, this number is 1 minus half of the P  value associated with −2 log likelihood difference. The result is shown in figure 12.

Figure 12shows the cumulative probability that the parameter is less than the value in the x-axis. Thus, the probability that the parameter is less than 0 (i.e. , one would live few days with the alternative chemotherapy) is 13%. Because 5 is the nadir, there is a 50% probability that the value of the parameter is less, and a 50% probability that the parameter is greater than 5. However, what if the patient wants to live 10 yr? There is an 87% chance that the true value of the parameter is less than 10 yr. There is only a 13% chance that the true value is greater than 10 yr, based on the data.

Standard reporting methods typically test a parameter to see if it is different than 0 at a P  value of less than 0.05. This method is woefully inadequate because the question of interest is really the following: “What is the likelihood of the parameter having a magnitude of interest?” Similarly, peer reviewers often ask for a power analysis, which is not particularly useful after the study has been completed. Most likely, the information they are really seeking is what we can infer about the possible values of the parameter from the study data. If the log likelihood profile is very flat, that says that we learn almost nothing about the parameter from the study data. If the log likelihood profile is quite narrow, then we know a great deal about the parameter. In either case, conversion of the log likelihood profile to the cumulative distribution allows visualization of the range of values of the parameter. The probability of the parameter being less than any clinically meaningful value along the x-axis is the cumulative distribution function (CDF) shown on the y-axis. Conversely, the probability of the parameter being above any clinically meaningful value is 1 − CDF.