An open-source framework for physiologically-based pharmacokinetic modeling of kinetics in the female reproductive tract.

Medical research and drug development have historically prioritized male physiology, leaving critical gaps in understanding female-specific responses. Adverse events in the 1960s and 1970s related to drugs such as thalidomide and diethylstilbestrol or devices like the Dalkon Shield ( Adam et al., 1977 ; Hoover et al., 2011 ; Mumford and Kessel, 1992 ) led to the exclusion of women from clinical trials as an attempt to protect women. Further, concerns about hormonal variability, and the longstanding assumption that data from men can be generalized to women resulted into a bias in medical research and drug development ( Merone et al., 2022 ; Holdcroft, 2007 ). Since the early 1990s, however, the inclusion of women in clinical trials has been generally required 1 , recognizing that drug pharmacokinetics (PK) and pharmacodynamics can differ between sexes ( Pare et al., 2011 ; Vasisht et al., 2021 ). Hence, adequate female representation is considered important to detect clinically relevant sex-related differences ( Merkatz and Junod, 1994 ; Gross et al., 2022 ). A lack of data stratified by biological sex can compromise drug effectiveness and increase the risk of adverse effects due to physiological variation ( Oi Yan Chan et al., 2022 ). Nevertheless, the hesitations to include women in clinical trials are still present ( Rehan et al., 2024 ). Model-informed drug development can support planning for safer and more targeted trials, making the inclusion of women in early phases both feasible and justified. Additionally, model-informed drug development supports sex-specific dosing strategies and, when combined with in vitro methods, helps fill data gaps in situations where direct studies may not be possible. Physiologically based pharmacokinetic (PBPK) models are a well-established mechanistic modeling framework to predict a compound's absorption, distribution, metabolism, and elimination. They allow the investigation of variability sources originating from physiological and anatomical differences, by population or individual-driven model parameterization 2 , 3 . They have also been applied to study sex-specific risk for chemical exposure ( Xiao et al., 2022 ; Hartmanshenn et al., 2016 ). Even though PBPK models are well-positioned to perform sex-specific analysis, very few modeling frameworks provide a sex-specific model structure, meaning most models do not represent the tissues of the reproductive organs or the breast. The female reproductive tract (FRT) consists of external organs (labia majora and minora, clitoris, vaginal opening, hymen and urethra) and internal organs (vagina, cervix, uterus, fallopian tubes, ovaries). It plays key roles in reproduction, hormonal function, and sexual health. The anatomical and physiological features of the female reproductive system give rise to a distinct set of healthcare needs — encompassing vulnerability to infections, conditions such as endometriosis and pelvic inflammatory disease, various cancers, hormone-related disorders, menstrual issues, fertility challenges, and menopausal complications. Representing the organs of the female reproductive tract ( Figure 1 ) explicitly in a PBPK model offers several benefits. Firstly, such a model extension enables the prediction of tissue and fluid concentrations within the FRT, after both local and systemic exposure. Secondly, it provides valuable predictions for applications where the FRT is the primary treatment target. Thirdly, establishing tools to support the development of vaginal and uterine drug administration systems—both of which bypass first-pass liver metabolism—can aid drug development and clinical trial design. However, despite the advantages of drug delivery via the FRT ( Sanchez Armengol et al., 2025 ), the development of such local administration systems is hindered by the lack of early drug discovery tools, the high cost of experimental set-ups, ethical concerns, and limited patient consent for clinical studies ( Adams et al., 2025 ). Representation of the Physiologically based Pharmacokinetic model structure. The left figure gives a schematic representation of the organs of the FRT. Furthermore, the different routes of administration that have been included in the exposure simulations are illustrated (Created in BioRender) 8 . The right figure gives a schematic overview of the full-body PBPK model. The presented extension module includes five tissue compartments (endometrium, myometrium, vagina, cervix and breasts), marked in pink, and two fluid compartments (uterine fluid and cervicovaginal fluid), presented in white. The PK-Sim base PBPK model already includes all the remaining organs. The literature around FRT PBPK models is limited. ( Thakur et al., 2024 ) and ( Kay et al., 2018 ) published PBPK models that comprehensively described the physiology and anatomy of the female reproductive tract organs. Both studies highlight the utility of these models for conducting cost-effective and accessible in silico mechanistic studies of vaginal drug absorption. Despite this potential, the studies also underscore that data gaps remain a major challenge when building such female-tailored PBPK model extensions. In this work, we developed a PBPK extension module for the female reproductive organs available in the Open Systems Pharmacology (OPS) Suite 4 and it is available on GitHub 5 . This module can extend any PBPK model available in the OPS Suite by adding the FRT organs. The mechanistic FRT extension module enables the prediction of tissue and plasma concentration following oral and intravenous administration as well as vaginal and intrauterine exposure. It was verified by comparing model predictions to clinical data for levonorgestrel and metronidazole, considering in total four routes of administration (intravenous, oral, vaginal and intrauterine delivery). Predicting drug exposure in female reproductive tract tissues is important for applications such as contraception, reproductive toxicology, and the treatment and prevention of communicable diseases. Hence, a robust extension module as part of a larger modeling ecosystem is expected to enhance the applicability of PBPK models to account for female-specific needs during drug development, clinical trial designs, and risk assessment.

A Shiny application (available at https://femsim.esqlabs.com/ ) was developed as a user-friendly and interactive tool for exploring the Female Reproductive Tract PBPK model, eliminating the need for prior experience with PBPK software. The app allows users to simulate local FRT concentrations for levonorgestrel and metronidazole, adjust dosing scenarios, and visualize predicted exposure in key reproductive tissues and compartments. The user can see the differences of simulated concentrations on the selected tissues over time for different administrations within all possible combinations of tissues and selected administrations for each drug. The FemSim Shiny App, is dependent on the {esqlabsR} package, for simulating the different scenarios selected by the user. It records no data from the user, and all the calculations are done on the fly (within the RAM). Once the calculation is complete, the {dygraphs} package 12 is used to plot the concentrations on different tissues over time within an interactive plot, where users can select different tissues after different drug administrations, change the Y-axis units and zoom in to any time of interest.

PBPK models were developed using PK-Sim® and MoBi® version 12.0 6 . Initial compound models were created in PK-Sim and extended in MoBi to incorporate the developed FRT extension module. Parameter identification was performed in MoBi. Data extraction from the published literature was done using WebPlotDigitizer v4.8. Sensitivity analysis and visualization of simulation results were performed in Claude (Anthropic) was used to assist in the preparation of R code for figure formatting and visualization. Simulations were performed using the CVODE solver with standard solver setting: absolute tolerance of 1E-10 and relative tolerance of 1E-5. The {esqlabsR} package 7 was used to run simulations, calculate pharmacokinetic parameters and generate graphical outputs. The PBPK model builds up in the standard whole-body PBPK model structure in PK-Sim ( Figure 1 ). The standard human model comprises 18 organs and tissues. Each organ is represented as a separate compartment, which is divided into the following four sub-compartments: plasma, blood cells, interstitial, and intracellular space ( Willmann et al., 2003 ). All absorption, distribution, metabolism and excretion processes are captured by ordinary differential equations ( Kuepfer et al., 2016 ). The PK-Sim parameter database enables model parameterization adjusted to individual characteristics such as age and biological sex. To represent the female anatomy a breast compartment was added during the model extension process. The FRT structure is developed as a separate extension module that can be added to any PBPK base model in MoBi version 12.0. The model extension ( Figure 1 ) simplifies the complex anatomy of the FRT, since there was no sufficient data to characterize all tissues of the FRT. The module consists of four FRT tissue compartments (endometrium, myometrium, vagina, and cervix) and two fluid compartments, uterine fluid and cervicovaginal fluid. A breast compartment is also included in the module to improve the overall representation of female physiology, although the focus of this work remains on the organs of the FRT. All tissue compartments have the same sub-compartmentation as the tissues of the base whole-body structure. Each tissue is connected via the systemic circulation, allowing drug exchange between tissues. Direct transfer between the endometrium and myometrium, despite their anatomical adjacency, is neglected due to the lack of data to quantify interlayer permeability. The whole-body PBPK base model parameters are obtained from the PK-Sim database that provides age and sex-dependent model parameters 9 . A structured literature search was conducted to identify anatomical and physiological data relevant to the female reproductive tract. Studies were included if they reported quantitative measurements in humans, focusing on adult, non-pregnant women. Data from animal studies was not included due to unclear translatability. The derived parameters include compartment volumes and specific blood flow rates, as presented in Table 1 . Data for detailed characterization of tissues compositions for the FRT tissues was not available. Therefore, we adopted the parameterization of the muscle compartment for those tissues. The tissue-specific parameterization of breast compartment, referring to the total breast tissue, was adopted from the pregnant population in the OSP Suite database ( Dallmann et al., 2017 ). Overview of physiological parameters for the female-specific organs of the extension module. Assumption based on variability observed for other tissues. Compounds released into the uterine or cervicovaginal fluid enter systemic circulation by crossing the tissues of the FRT. The transport equations for the compound transfer between the uterine fluid and the endometrial tissue, as well as between the cervicovaginal fluid and vaginal and cervical tissue are based on the mathematical description of a diffusion process. The equation describing the amount of compound in the vaginal or uterine fluid N f l u i d reads: d N f l u i d d t = D   S A t i s s u e H t i s s u e   C t i s s u e · f u t i s s u e · s c a l i n g − C f l u i d · f u n i o n + d N f o r m d t . (1) D is the diffusion coefficient in the tissue through which the compound diffuses. Since experimental data on drug transport across tissues of the FRT organs were unavailable, we estimated the diffusion coefficients for the cervix and vagina using a Quantitative Structure-Activity Relationship (QSAR) model developed by ( Chen et al., 2015 ). Although designed initially for skin tissues, we followed the rationale of Thakur et al. assuming its applicability to vaginal tissue due to physiological similarities. The diffusion coefficient for the uptake from the uterine fluid was initially assumed to be equal to the value for cervical and vaginal uptake. However, numerical optimization was needed to ensure that the model captured the clinical data appropriately. The dimensions of the tissues are described by a surface area (SA) and a thickness (H). C denotes the concentration in the tissue and fluid, respectively. The unbound fraction in tissue is denoted by f u t i s s u e . As only unionized compounds can diffuse into the tissues, this fraction is captured by f u n i o n . The fraction unionized was derived using the Henderson–Hasselbalch equation ( Po and Senozan, 2001 ), which relates the drug’s pKa to the local pH and thus allows calculation of the proportion of molecules in the unionized form. We also included a scaling factor in the transport formula. It can be interpreted as a correction factor for the unbound fraction in tissue and can be optimized when tissue concentration data is available. By default, this value is set to one but provides one degree of freedom for numerical optimization. In practice, the unbound fraction in tissue is rarely reported. Hence, the parameter is usually derived based on assumptions. The introduced scaling factor can be interpreted as a correction factor for additional and unknown binding in the tissues, which has not been captured when deriving f u t i s s u e . In case of local administration, e.g., an intrauterine device or a vaginal tablet, the source term d N f o r m d t describes the mass transfer from the formulation into the uterine or cervicovaginal fluid. The absorption model presented in Equation 1 assumes that all compartments are well-stirred, and that the fluid-tissue interface is uniform. While d N f o r m d t is zero when working with systemic exposure routes, the compound release d N f o r m d t from an intrauterine or intravaginal formulation to the fluids of the FRT can be formulated with a release rate or a release profile. The most appropriate mathematical expression for the release term depends on the administration route, formulation, and available data. For this work, we evaluated two expressions for d N f o r m d t , namely, a zero-order process with a release rate k r e l e a s e (mg/h): d N f o r m d t = k r e l e a s e (2) and a Weibull function: d N f o r m d t = 1 − exp − t b a ⁡ . (3) The parameters a and b characterize the shape and the scale of the Weibull release profile. For the levonorgestrel IUD 52 mg longitudinal release data as well as an expected average release rate along with PK data ( Cicali et al., 2021 ) are publicly available. Therefore, we investigated the differences in predicted plasma concentration when using (A) a zero-order release term ( Equation 2 ) and (B) a Weibull release term ( Equation 3 ) for the Levonorgestrel (LNG) IUD 52 mg. The compounds were selected based on the availability of in vivo plasma, tissue and fluid concentrations after systemic and local administration to allow the development and evaluation of a base PBPK model. Furthermore, clinical data for the characterization of the different local drug administration formulations needed to be available. The selected compounds were levonorgestrel and metronidazole. Levonorgestrel is a synthetic progestogen widely used in various hormonal contraceptive formulations. Cicali et al. (2021) developed a PBPK model for orally administered levonorgestrel, available on GitHub. Their model accounts for CYP3A4-mediated metabolism as well as binding to systemic albumin and sex hormone-binding globulin. Their drug-specific PBPK model parameters are summarized in Table 1 in the Supplementary Material . We used the PBPK model from Cicali et al. as a base model and extended it with our FRT module to predict levonorgestrel concentrations in female reproductive tract tissues following both oral and intrauterine administrations. The simulated dosing strategies were set up to match the study protocols of Cicali et al., Nilsson et al. and Apter et al. Compound-specific PBPK model parameters that were introduced for the FRT extension are provided in Table 2 . Thakur et al. (2024) used a QSAR established for skin ( Chen et al., 2015 ) to predict the diffusion coefficient for the vaginal tissue. We extended the application of the QSAR to have an initial parameter estimated for the cervical and uterine tissue, knowing that their rational for its applicability did not hold for uterine tissue. Therefore, to predict the tissue concentration reported by Nilsson et al., the diffusion parameter for the uptake of the uterine tissue was refitted. Furthermore, we calculated the scaling factor of the unbound fraction in endometrial and myometrial tissue. Our interpretation of this scaling factor is that it should account for differences in LNG target availability, as reported for these tissues ( Thi et al., 1975 ). For the intrauterine device containing 63 mg (IUD 1) and 52 mg (IUD 2) LNG, we characterized release functions. A Weibull function could not be characterized for IUD one due to the lack of longitudinal release data. The parameters of the corresponding release terms are also given in Table 2 . Overview of levonorgestrel-specific parameters for the FRT module. Results of a multiple start optimization are presented in Table 6 and Figure 3 in the Supplementary Material . Metronidazole is an antibiotic and antiprotozoal medication used to treat a range of infections, including bacterial vaginosis, trichomoniasis and endometritis. A PBPK model for IV administration of metronidazole was taken from the literature ( Nil et al., 1982 ). The calculation methods were Rodgers and Rowlands for partition coefficients and PK-Sim Standard for cellular permeabilities. The model was transferred to PK-Sim/MoBi V12 and extended for oral and vaginal administration. The drug-specific PBPK base model parameters as derived by Dallmann et al. For the oral formulation, a dissolved formulation was applied to the immediate release formulation, while a Weibull was fitted to clinical data for the extended-release formulation. The parameters specifically for oral demonstrations are in Table 3 in the Supplementary Material . Metronidazole-specific parameters for different administration types. All release profiles are described with Weibull functions. Levonorgestrel and metronidazole were selected as validation compounds because of the availability of plasma and tissue PK data for multiple systemic and local routes of administration. The primary objective of the validation was to evaluate model performance across different administration routes rather than across physicochemical space. It should be noted that generalizability to compounds with markedly different physicochemical properties, such as weak bases with pronounced pH-dependent ionization or highly lipophilic molecules, remains to be established. To validate the model’s predictions, we collected plasma and tissue levonorgestrel concentrations following both oral and intrauterine administrations. To replicate the results of Cicali et al. (2021) , we included their dataset, which reports levonorgestrel plasma concentrations after a 0.09 mg daily oral dose. We also incorporated the study by Nil et al. (1982) , which details plasma and tissue (endometrium, myometrium, fat, and breast) concentrations after 7 days of 0.25 mg daily oral administration and 36–49 days of intrauterine administration with a device releasing 30 μg of levonorgestrel per day. Finally, the ( Apter et al., 2014 ) study on levonorgestrel plasma concentrations for various IUD formulations over extended period of time allowed us to explore different mathematical models for drug release from the IUD and their impact on plasma concentration predictions. For metronidazole, clinical data were collected from the literature for intravenous, oral and vaginal administration. The PBPK model from Dallmann et al., was first developed and compared to PK data of healthy women after 500 mg IV administration ( Houghton et al., 1979 ), and further evaluated in pregnant women receiving the same dosing regimen ( Visser and Hundt, 1984 ). Model performance was assessed by comparing predicted plasma concentration profiles to observed clinical data ( Shepard et al., 2015 ; Kuemmel et al., 2020 ). Visual comparisons were made by overlaying clinical data with the predictions. The predicted-to-observed ratios were computed, and 2-fold error interval, as commonly applied ( van Borselen et al., 2024 ), was used to assess model acceptability. Furthermore, PK parameters such as C max and AUC 0-t were calculated for the model predictions and compared to observed values from the study reports. The (GMFE) G M F E = 10 ∑ log ⁡ P p r e d P o b s n (4) was calculated as a quantitative measure to assess model performance based on predicted and observed PK parameters, P p r e d and P o b s respectively. For data set with single point measurements, it is not feasible to calculate C max and AUC 0-t . In such cases, the observed and predicted concentration, C o b s and C p r e d respectively, are used to derive the prediction error (PE) in percent P E = C p r e d − C o b s C o b s   · 100 % . (5) A local sensitivity analysis was performed to assess the impact of model input parameters on the model’s predictions and build confidence in our modelling approach. Furthermore, the results can identify where experimental data would most effectively support the modeling process. The analyzed parameters included organ-specific blood flow and organ composition to understand better tissue-specific parameters introduced for the female reproductive tract organs. Parameters essential for vaginal and uterine absorption models were also evaluated. In total, 40 parameters have been analyzed, and a complete parameter list is provided in Table 5 in the Supplementary Material . The effects on the predicted time curves as well as PK parameters (AUC, C max , t max ) for plasma and tissue (breast, vaginal, cervical, uterine, myometrium and endometrium) concentration were considered in the sensitivity analysis. The analysis was carried out using the default implementation of a local Sensitivity analysis in {esqlabsR} package. Exploratory population simulations were performed to assess the impact of variability in the female population. Female European populations (n = 100) were created to match study demographics for age, weight and height. If no data was available, a population was created with age ranging between 29 and 31 years, incorporating the default population variability in OSP. Furthermore, additional variability was implemented for the FRT organs’ volume and specific blood flow. Figure 1 shows the mean and standard deviation for organ volume and specific blood flow collected from literature. A lognormal distribution was assumed for all FRT parameters in the population.

The established PBPK models were applied to predict levonorgestrel and metronidazole concentrations in plasma and tissues following both systemic and local administrations. Additionally, we evaluated how the selection of the drug release function influences exposure predictions. Lastly, a local sensitivity analysis was conducted to explore the impact of tissue-specific parameters and parameters of the vaginal/uterine absorption model on model outcomes. The parameters of the base PBPK model were not refined during model extension. We, however, performed model fitting for the parameters of the uterine absorption process. Table 2 provides an overview of the absorption process parameters presenting their initial and optimized values. The uterine diffusion coefficient was refitted using concentration data from endometrium and myometrium. A bracketing uncertainty analysis for the uterine and oral LNG administration (±1 log 10 around the QSAR-derived estimate) was performed to assess the impact of QSAR transferability uncertainty on model predictions ( Figures 4 , 5 in the Supplementary Material ). For oral administration, the diffusion coefficients have no significant impact on exposure predictions, as systemic exposure is driven by blood flow. For intrauterine administration, the bracketing analysis revealed asymmetric uncertainty for the uterine diffusion coefficient. The upper bound had a modest impact on predicted tissue concentrations, while the lower bound resulted in substantially reduced predictions. The physiological parameters of the base PBPK model were those of a 30-year-old healthy female individual as characterized in the PK-Sim database. The mean parameters for the organs of the female reproductive tract are given in Table 1 . The predicted plasma and FRT tissue concentrations following oral ( Figure 2A ) and IUD ( Figure 2B ) administration are presented in Figure 2 . The predicted vs. observed plots ( Figures 2C,D ) visually indicate that the model fits the data well. To quantify the model prediction capabilities, the prediction error (PE) ( Equation 5 ) was calculated for the eight observation-prediction pairs ( Table 4 ). Six out of eight PE were under 50%, indicating a very good prediction and two had a PE under 100% marking a good prediction. The direction of the PE indicated that the model tends to slightly underpredict. Due to the sparsity of the data, it was not feasible to calculate the AUC or C max. Hence, the GMFE cannot be provided. Simulated concentration–time profiles following (A) daily oral administration of levonorgestrel (250 μg) and (B) intrauterine device administration with a daily release rate of 30 μg. Observed plasma and tissue concentrations reported by Nilson et al. are shown as symbols. Solid and dashed lines represent model-predicted concentrations, with colors corresponding to the respective observed data. Observed versus predicted concentrations for the oral formulation are shown in (C) , and corresponding results for the IUD are presented in (D) . Summary of observed and corresponding predicted concentrations for plasma, fat tissue, endometrium and myometrium after oral and IUD administration of LNG, along with their prediction error. To capture the tissue concentration in both the endometrium and myometrium, the scaling factor introduced in Equation 1 was optimized. The fitted values were 0.11 for the endometrium and 0.23 for the myometrium, resulting in reduced availability of freely moving levonorgestrel in both tissues. The most significant reduction can be observed in the endometrium. This is qualitatively in line with experimental observations, which suggest that progesterone receptor abundance is higher in the endometrium compared to the myometrium ( Thi et al., 1975 ). Furthermore, these findings align with our hypothesis that the scaling factor accounts for differences in target availability between tissues. The release from a drug administration device, such as an intrauterine device, can be described mathematically in different ways. The choice of release model depends primarily on the available device-specific release data and clinical data for model calibration. For the IUD 52 mg LNG, we compared two release models: (A) a zero-order release and (B) a Weibull function. The Weibull parameters ( Table 2 ) were obtained by fitting the function to reported release data (see Supplementary Figure 2 ). In both cases, we used the release parameters derived from the reported release experiments without further refinement. The objective was not to identify an optimal release model, but to evaluate how different mathematical descriptions of device release influence simulated plasma concentration, and to identify where current approaches perform adequately and where refinement would be needed in future work. Figure 3 shows the predicted plasma concentration–time profiles for the two release models. Visual inspection suggests that the Weibull function captures the overall shape of the release profile better than the zero-order model, particularly the initial rise and the gradual decline in concentrations over time ( Figure 3B ). However, pharmacokinetic parameters ( Table 5 ) indicate that the zero-order model better predicts overall drug exposure (AUC). The Weibull model, parameterized from fitting to estimated in vivo release rates (Mirena Labeling), systematically underpredicts plasma concentrations resulting in an underestimated AUC. Refining the Weibull parameters could improve this discrepancy. Neither release model predicted t max accurately; the zero-order model approximated C max but at an incorrect time point. Simulated plasma concentration–time profiles of levonorgestrel (LNG) following intrauterine device (IUD) administration, comparing two different release models: (A) zero-order release, assuming a constant daily release rate, and (B) Weibull release function, representing a time-dependent, decreasing release rate. Summary of observed and predicted LNG plasma concentration after IUD administration. The release from the device was described by (A) a zero-order release term and (B) a Weibull function. The impact of the release term on the prediction accuracy is captured by corresponding GMFE. In conclusion, the Weibull function leads to a better description of the qualitative feature of the release profile. However, in scenario where one is interested in exposure over longer periods, the zero order release term may be sufficient. Figure 4 shows the concentration-time profiles for metronidazole after intravenous, oral and vaginal administration. The predicted versus observed plots are included in the Supplementary Material . The GMFE on the Cmax in plasma was 1.39, 1.11 and 1.11 for intravenous (n = 3), oral (n = 6) and vaginal (n = 5) administration respectively. The GMFE on the AUC in plasma was 1.18, 1.10 and 1.14 for intravenous, oral and vaginal administration respectively. The prediction error for metronidazole concentration in myometrium and endometrium were −6% and −12% for intravenous and −24% and −27% for oral administration, respectively. PK parameters (AUC and Cmax) are listed in Table 6 , including predicted/observed ratio. Simulated concentration–time profiles for metronidazole following (A) a single intravenous administration of 500 mg (B) 400 mg thrice daily oral administration (C) single vaginal gel administration of 37.5 mg and (D) 400 mg twice daily oral administration. Mean observed plasma and tissue concentrations are shown as symbols, with error bars representing the standard deviation. Solid and dashed lines represent simulated concentrations, with colors corresponding to the respective plasma or tissue concentrations. Summary of calculated PK parameters for metronidazole. Using local sensitivity analysis, we investigated the model parameters’ impact on AUC, C max and t max after uterine LNG and vaginal MET administration. Figure 5 shows an overview of the most sensitive parameters for AUC following local drug administration levonorgestrel (A) and metronidazole (B). Results of the local sensitivity analysis for AUC. (A) Shows the results for the levonorgestrel model following uterine administration and (B) For metronidazole following vaginal administration. Only parameters leading to at least 1% change in AUC are presented. The investigated parameters were scaled with a factor in the range [0.1–10.0]. Color intensity indicates the magnitude of AUC change in percent. For vaginal administration, the most sensitive parameters relate to the vaginal and cervix tissues, while for uterine administration endometrial and myometrial parameters have a noticeable impact on model predictions. The analysis of the other PK parameters leads to the same conclusion. A representative population simulation for a female European population (n = 100) is shown in Figure 6 . The simulation was performed for a population with age ranging between 22 and 31 years, and weight varying between 51 and 76 kg according to the study demographics. Population simulation. A population simulation (n = 100) was performed for vaginal application of metronidazole 500 mg twice daily vaginal. The line represents the median plasma concentration and shaded area represents 5th to 95th percentile. Points show the data from Salas-Herrera 1991b.

Predicting drug exposure in female reproductive tract tissues is important in drug development, with applications in contraception, reproductive toxicology, and the treatment and prevention of communicable diseases. The physiology of the FRT organs, such as tissue volume, blood flow, and local pH, can influence how a compound is absorbed and distributed, making a mechanistic modeling framework like PBPK well suited to describe tissue pharmacokinetics based on organ-specific characteristics. Here, we present a publicly available PBPK model extension applicable to all PBPK models in PK-Sim/MoBi version 12.0, evaluated with levonorgestrel and metronidazole across at least two routes of administration each. The model was able to capture the pharmacokinetics of both compounds across different routes of administration well. The following discussion addresses the key assumptions, limitations, and areas for future refinement. We parameterized the extension module using publicly available data, focusing on quantitative data originating from human studies, and published PBPK models including the FRT organs from other PBPK software ( Thakur et al., 2024 ; Kay et al., 2018 ). As Table 1 indicates by the number of references supporting model parameters, there are generally few sources available to parametrize the model compartments. Furthermore, several sources did not provide sufficient information to characterize parameter variability, which limits the quantification of parameter uncertainty and represents a drawback for population simulation. Lastly, adequate tissue composition data were not available to characterize the endometrium, vagina, and cervix, leading to the use of muscle as a template organ. This assumption directly impacts the calculated partition coefficients for these tissues. To assess the influence of this uncertainty on model output, partition coefficients for all three organs were included in the local sensitivity analysis, which confirmed that they can impact model predictions. Based on these results, expected variation in tissue composition will affect predicted exposure in the FRT compartments. This therefore represents a source of model uncertainty that could be reduced by closing this data gap. The model does not include the fallopian tubes or ovaries as explicit compartments, due to insufficient human data in the public domain to support parameterization. Given the lack of human data, one could explore the data availability from non-human primates and evaluate their applicability to parameterize a human model. The sensitivity analysis identified vaginal and cervical tissue parameters as most influential for vaginal administration, and endometrial and myometrial parameters for uterine administration. This is relevant in the context of prediction uncertainty: parameters identified as highly sensitive that are insufficiently supported by data (see Table 1 ) represent the primary source of model uncertainty, and targeted experimental data collection for these parameters would be the most effective way to reduce it. Specifically, the model would benefit from data to characterize the cervical blood flow as well as data that allows a better characterisation of the tissue properties of the different tissues of the FRT. The current model parameterization captures the characteristics of a healthy female of child-bearing age best, but it does not incorporate menstrual cycle–driven changes in tissue volume or blood flow. Given the results of our sensitivity analysis, including menstrual cycle-driven changes would add value for certain applications. For example, changes in the endometrium volume and the uterine blood flow could improve certainty for uterine drug delivery systems. Overall, the data used to derive the physiological parameters presented in Table 1 was quite variable and the model parameterization could be further refined with additional data. The fluid volumes and their pH are assumed to be constant. Hence, the pH gradient observed for the cervicovaginal fluid is neglected. For levonorgestrel and metronidazole, this simplification did not have a significant effect on the fraction of unionized drug available for uptake from the vaginal fluid. For compounds where the unbound fraction is more sensitive to the pH range of the vaginal tract, sub-compartmentalization should be considered. Assuming a constant fluid volume rather than a dynamic turnover description is a further simplification, as secretion and drainage processes influence PK. This limitation is most consequential for drugs with short therapeutic windows or dosing regimens sensitive to the timing of luminal clearance. Calibrating a turnover rate would require additional data such as secretion rate measurements and longitudinal cervical-vaginal fluid concentration profiles with sufficient resolution to separate physical washout from mucosal absorption. Nevertheless, given the modular model development concept, these parameters and mechanisms can be updated as data become available. A major challenge in model development was the parameterization of drug transfer from vaginal and uterine fluids to target tissues. Ideally, in vivo release data would be used to parameterize this process, but in the absence of such data, in vitro release data generated under biorelevant conditions could serve as an alternative. However, an established framework for in vitro–in vivo extrapolation (IVIVE) for vaginal and uterine formulations is currently lacking. A QSAR derived for skin tissues was used for the parameterization of diffusion coefficients for the FRT compartments. While this represents the best available approximation in the absence of a tissue-specific QSAR framework, the application is not ideal given the differences between the tissues. The uterine diffusion coefficient was refitted to available concentration data, which partially mitigates the QSAR limitation for that compartment. We adopted a mixed approach, combining data-driven model fitting with a priori knowledge. A key challenge was the incompleteness of publicly available datasets. Notably, the limitation in dataset reporting plasma, tissue, and fluid concentrations simultaneously hindered fully data-driven optimization and has direct implications for parameter identifiability. In the absorption equation, the diffusion coefficient and the scaling factor serve structurally distinct roles, D governs the overall magnitude of flux symmetrically, whereas the scaling factor applies exclusively to the tissue-side driving force, and surface area and membrane thickness were fixed based on literature morphometrics. Nevertheless, fitting both D and the scaling factor to limited concentration data introduces a source of non-identifiability, and alternative parameter sets may describe the data equally well. The fitted parameters should therefore be interpreted as consistent with the available data rather than as uniquely determined mechanistic values. Alongside quantitative data, we leveraged qualitative insights to inform and constrain the model. For example, the scaling factor introduced to account for LNG availability in the endometrium and myometrium was supported by animal studies showing higher progesterone receptor prevalence in the endometrium. Our fitted value aligned with this, indicating a smaller scaling factor for the endometrium, consistent with a larger unbound fraction of LNG and potentially greater receptor binding in that tissue. In conclusion, despite data limitations and the inherent constraints on parameter identifiability, the simulation results indicate that this mixed approach was an appropriate and pragmatic strategy for parameterizing the model extension. Overall, closing identified data gaps for tissue characteristics, populations information and local transport information would substantially reduce this uncertainty in future work. Similarly, data to parameterize the vaginal formulation of metronidazole were not available from the literature. Therefore, various vaginal applications of metronidazole were modelled assuming zero-order release kinetics, and automated parameter identification was applied. Furthermore, the administered dose was scaled to align the PBPK model predictions with observed clinical bioavailability. The reduced apparent exposure can be partially explained by the self-cleansing action of the vaginal environment, which promotes drug removal from the administration site. Nonetheless, other less-understood physiological (e.g., local metabolism in the female reproductive tract) or formulation-related effects likely play a role, indicating the need for a better understanding of the female reproductive tract. For the intravenous scenarios, some larger deviations in Cmax were observed, which may reflect uncertainties related to sampling design, formulation assumptions, and tissue partitioning rather than fundamental model misspecification. The Shiny app offers a practical way for researchers, students, and developers to engage with the model immediately and test different assumptions or use cases. By lowering technical barriers and enabling rapid experimentation, it complements the webinar and online training modules and helps broaden the community of users able to work with the reproductive tract model. It also supports early hypothesis generation, teaching, and potential research extensions, reinforcing the long-term value of an openly accessible modelling framework. To conclude, a mechanistically informed extension module for the female reproductive tract was developed and implemented within the Open Systems Pharmacology PBPK framework. Integration of this module into existing PBPK models for levonorgestrel and metronidazole enables descriptions of pharmacokinetics following both local and systemic administration. The module is openly available to the scientific community via GitHub, and a user-friendly Shiny app allows exploration of plasma and tissue concentrations for these drugs with minimal technical requirements. Overall, the presented work provides a valuable resource to address female-specific considerations in drug development, clinical trial design, and risk assessment. Valuable model extensions for future work include broadening the chemical space the model is applied to, incorporating mechanistic detail such as dynamic fluid turnover, drug-protein binding and pH stratification in the vaginal tract, and accounting for cycle-driven variability and population-specific parameterizations, for example, in pregnant and menopausal women.