Breastmilk is an essential and valuable source of sustenance and immunological protection for infants, providing important nutrients and antibodies that promote strong growth and development (Wedekind and Shenker, 2021; Czosnykowska-Łukacka et al., 2018). Over 50% of women take medication during the breastfeeding period, but data regarding the potential exposure of their nursing infants to maternal medication are typically limited, forcing women to choose between their own health and potential harm to their children (Bérard and Sheehy, 2014; Stultz et al., 2007; Al-Sawalha et al., 2016; Ojara et al., 2023; Wang et al., 2017). Although there is increasing demand and regulatory recommendation of the inclusion of lactating women in clinical trials, there is still much work to be done to support safe medicine use during breastfeeding (FDA, 2019; EMA, 2024).

The inclusion of breastfeeding women in clinical trials does present with ethical issues, and so the development of alternative methods of to support increased understanding of potential infant exposure is desirable (Nauwelaerts et al., 2021). A number of non-clinical methods to assess the passage of maternal medication into the breastmilk are available, including in vitro human (Kimura et al., 2006; Andersson et al., 2017) and animal models (Ito et al., 2015; Al-Bataineh et al., 2009), in vivo animal models (McNamara et al., 1992; Osz et al., 2017) and in silico models for the prediction of milk concentrations (Begg and Atkinson, 1993; Fleishaker et al., 1987; Koshimichi et al., 2011; Nauwelaerts et al., 2023). Differences in transporter expression of cell lines have been noted within the literature and in vitro work is also complicated by the use of cell-specific media over natural biological matrices (Nauwelaerts et al., 2021; Qu et al., 2015). Similarly, species-specific differences may limit the use of animal models. Rodent models, in particular, have shown a much higher milk-to-plasma (M/P) ratio for some drugs than has been reported in humans (Ito et al., 2013). More recently, Gottingen Minipigs have shown promise as an in vivo model, with a strong correlation between minipig and human M/P ratios being reported [unpublished data from IMI ConcePTION (Annaert, 2024)]. Full reviews of in vitro and in vivo models have been published previously (Nauwelaerts et al., 2021; Ventrella et al., 2019).

The use of mathematical modelling for the in silico prediction of maternal M/P ratios is of interest as it removes the concerns regarding species-specific differences in lactation mechanisms. Such models have been used to successfully predict the M/P ratio for medications including primaquine, theophylline, ondansetron and sotalol (Begg and Atkinson, 1993; Pan et al., 2023; Abduljalil et al., 2022; Job et al., 2022). Quantification of maternal plasma concentrations through clinical study and/or physiologically based pharmacokinetic (PBPK) modelling allows M/P ratios to be translated into milk concentrations and therefore the subsequent prediction of infant exposure to maternal medications through breastfeeding (Abduljalil et al., 2022; Job et al., 2022; Abduljalil et al., 2021; Pansari et al., 2022).

Mathematical models for predicting the drug transfer into milk have been of interest for decades, with Atkinson and Begg, and Fleishaker et al separately publishing mathematically identical models over 30 years ago (Begg and Atkinson, 1993; Fleishaker et al., 1987; Atkinson and Begg, 1990). These models quantify a ‘phase distribution’ theory of drug partitioning in which the unbound and unionised fraction of drug in the plasma and the breastmilk exists in equilibrium, and differences in milk and plasma concentrations arise due to differences in protein binding, ionisation and lipid partitioning between the two matrices. Natural variability in milk pH, fat content (creamatocrit) and protein content/binding therefore lead to prediction of variable M/P ratios, assuming these factors are sufficiently considered. This Phase Distribution model is supported by in vitro bioanalysis conducted by Atkinson and Begg, in which formulae for the prediction of milk protein binding [fu _m (Atkinson and Begg, 1988a)] and lipid partitioning [Papp _milk (Atkinson and Begg, 1988b)] at pH 7.2 were developed based on the known plasma protein binding (fu _p ) and octanol-to-water distribution coefficient (LogD) respectively. This formula for prediction of Papp _milk shows a steep relationship between LogD_7.2 and Papp _milk , but was developed using drugs with a LogD_7.2 of <3 only (See Supplementary Figure S2). Following this initial work, Atkinson and Begg employed linear regression to optimise predictions of acidic and basic drugs, with this ‘Log-Transformed’ model providing better predictions than the initial Phase Distribution model (Begg and Atkinson, 1993; Atkinson and Begg, 1990).

More recently, Abduljalil et al. incorporated these mathematical models into PBPK modelling software using an alternative formulation of the Phase Distribution model in which the lipid partitioning is represented by LogP_o:w (LogP_o:w model). This results in a significantly reduced predicted M/P ratio for lipophilic drugs but was shown to accurately predict the milk distribution of acetaminophen, alprazolam, caffeine and digoxin (Abduljalil et al., 2021; Zhang et al., 2022). Subsequent work by the same group has successfully used both the Phase Distribution and Log-Transformed model to predict infant exposure to a number of medicines (Pan et al., 2023; Abduljalil et al., 2022; Pansari et al., 2022). Although this LogP_o:w model represents a misinterpretation of the original Phase Distribution model, it is a useful comparator to study the way in which drug distribution into the milk lipid is incorporated into the model.

In best practice, use of these models includes a sensitivity analysis which accounts for the physiological variability of breastmilk (Pan et al., 2023; Pansari et al., 2022). This is particularly important because lactation studies do not typically include paired recording of factors which are predicted to affect the M/P ratio, such as pH and creamatocrit (Crt), which itself present challenges for model development. Despite their increasing use, to the best of our knowledge there has been no widescale review of the reliability of the Phase Distribution, Log-Transformed or LogP_o:w model accounting for the biological heterogeneity of human breastmilk. Additionally, the formulae for prediction of lipid partitioning and milk protein binding have come under little scrutiny since their conception. As such, the present paper aims to assess the reliability of existing mathematical models for predicting the M/P ratio of a diverse list of medications, and to subsequently interrogate additional factors that may alter prediction reliability.

All models evaluated showed misprediction of a significant proportion of the drugs. Of the pre-existing models investigated, the Log-Transformed model appeared to be most reliable for the prediction of drug distribution into milk for basic drugs, however the prediction of acidic drugs was very variable and there is no formula for predicting neutral medications. We subsequently investigated whether interrogation of the formulae underpinning the Phase Distribution model could improve predictability. An alternative model for prediction of milk lipid distribution (Bartels-LogD model) improved predictions, but the use of observed fu _m or milk pH values; alternative equations for predicting fu _m ; or incorporation of ionisation into milk protein binding predictions failed to significantly improve model reliability. OCT1 substrates and drugs with a low fu _p were identified as more likely to be underpredicted.

The limited predictability of the Log-Transformed model has been highlighted previously, where analysis of 69 medications, assuming a milk pH of 7.2 and a Crt of 4.5%, showed no correlation between the predicted and observed M/P values [R² = 0.01; (Larsen et al., 2003))] and this concern has also been raised by other groups (Yang et al., 2022; Pansari et al., 2024). Although the LogP_o:w model showed reasonable predictability compared with the other models, the use of LogP rather than the partition coefficient (i.e., the exponent of LogP; see Equations 4, 5) results in the model failing to recapitulate how changes in Crt affect drug distribution. For example, the LogP_o:w model predicts the M/P ratio of the lipophilic compound labetalol to increase from 0.80 to 0.91 as Crt increases from 3% to 12%, whereas the Phase Distribution model predicts this to rise from 2.95 to 9.50. As such, focusing on the non-log-transformed Phase Distribution model to derive better predictions was deemed most appropriate. This is supported by the knowledge that the LogP_o:w model is in fact a misrepresentation of the original Phase Distribution model [as detailed in (Zhang et al., 2022)], and that the authors have subsequently employed the Phase Distribution and Log-Transformed models (Pan et al., 2023; Abduljalil et al., 2022; Pansari et al., 2022). The ABI and Bartels-LogD models predict a significantly lower Papp _milk for drugs with a LogD of more than 3 (see Supplementary Figure S2), and the improved predictions of lipophilic drugs using these models demonstrate that the linear extrapolation of the Atkinson/Begg formula included in the Phase Distribution model may not be appropriate.

Our data showed that OCT1 substrates were significantly underpredicted, and this is consistent with a role for OCT1 in the active transport of drugs into the breastmilk. To the best of our knowledge, there has been little study into this role of OCT1 during lactation in humans, but rodent models have shown a reduction in the OCT1 substrate, thiamin, in the milk of OCT1 knockout mice (Kato et al., 2015). In support of this, OCT1 is upregulated during lactation (Kimura et al., 2006; Alcorn et al., 2002). Similarly, BCRP has been shown to be upregulated during lactation (Ahmadzai et al., 2022; Sychterz et al., 2024), and it was therefore surprising to see that the prediction reliability for BCRP substrates was not significantly different to that of drugs not mediated by transporters (Figure 4A; p = 0.068). We acknowledge that this may have been the result of using a limited dataset only, with other groups previously highlighting that BCRP substrates are likely to be actively transported into the breast milk (Ito et al., 2015; Yang et al., 2022; Sychterz et al., 2024; Gong et al., 2024). Similar to the present study, these previous works typically rely on a small number of BCRP substrates only. Prior work has demonstrated that in vitro to in vivo extrapolation (IVIVE) can be used to better predict the M/P ratio of BCRP substrates (Ito et al., 2015; Yang et al., 2022). This was considered outside the scope of the current work but is certainly an area in which more work should be focused, and gathering data from more substrates of BCRP (and other transporters) would clearly be beneficial. At present, these data are challenging to interpret given that the majority of medications are substrates for multiple transporters that are expressed within the breast. Our data agree well with the previously suggestion that MDR1 and MRP2 do not play a significant role in the transfer of drug into the breast milk (Alcorn et al., 2002; Ahmadzai et al., 2022; Gong et al., 2024).

From a clinical perspective, the underprediction of M/P ratios may contribute towards dosing nursing infants with unsafe quantities of maternal medications via the breastmilk, and therefore understanding the role of transporters such as BCRP and OCT1 is imperative. When considering only drugs that are not OCT1/BCRP substrates, 97% of medications were predicted within 2-fold of the observed data or overpredicted (2/64 underpredicted) when using the pH range of 6.63–7.62. Although overpredictions are not desirable per se, these scenarios allow the use of a model-based approach for estimating worst-case scenarios in the presence of limited clinical data. Whilst work should be encouraged to continue improving the reliability of such models, understanding of current limitations to ensure models are employed in appropriate ways is also beneficial.

It was somewhat surprising to see that using observed fu _m values and alternative models for predicting fu _m did not improve or even significantly affect predictions (Figure 3). The model assumes that only the unbound fraction of a drug is able to partition into milk lipids, and so increasing the milk protein binding will lead to a ‘compensatory’ increase in lipid partitioning, and therefore little change in predicted M/P ratio. Although this physiological mechanism supports our findings, there is little-to-no experimental evidence to verify the interplay between protein binding, lipid partitioning and ionisation. Atkinson and Begg’s original protein binding (Atkinson and Begg, 1988a) and lipid partitioning (Atkinson and Begg, 1988b) experiments were all performed at pH 7.2, as were more recent, similar, experiments (Ito et al., 2013). More detailed in vitro work is therefore needed to quantify the relationship between these factors in order to develop more reliable models.

The limited availability of paired milk pH, Crt, protein binding and M/P ratio data provides additional complications. The models can be used to predict the effect of variable milk physiology on drug accumulation, and therefore infant exposure, however there are little data available to verify this predicted effect of Crt or pH. Our publication list contained only two studies with an AUC M/P ratio and Crt data (Rampono et al., 2000; Kristensen et al., 2007). In one of these, fore- and hindmilk were analysed separately and shown to have Crt values of 6.2% and 13.7% respectively, which correlated with a 2.3-fold increase in mirtazapine concentration [LogP 2.9; (Kristensen et al., 2007)]. This change is recapitulated by both the Phase Distribution model (2.1-fold increase) and the Bartels-LogD model (1.89-fold increase), with the Bartels-LogD model more accurately predicting the M/P ratio (P/O ratio of 0.87–1.65 vs. 3.10–6.56; all at pH 7.12). Such a difference in the Crt of fore- and hindmilk is well established throughout the literature, with hindmilk typically exhibiting Crt values 2-3 times higher than that of the foremilk (Meier et al., 2006; Kristensen et al., 2007; Daly et al., 1993; Bowornkitiwong et al., 2023; Mizuno et al., 2009). These data may be beneficial when considering study design: if a drug has a low LogD at breast milk pH then it can be anticipated that samples of foremilk only will reflect drug concentrations in the whole milk. In contrast, for lipophilic drugs, sampling of whole milk or fore- and hindmilk will be considered more important. When we have a robust understanding of how changes in Crt affect dug distribution, it is plausible that modelling may be used to extrapolate whole milk values from fore-milk only, which may reduce the burden on women recruited to clinical trials. As stated, additional observed data with paired milk pH and Crt information are required to verify this suggestion, and will thus support the use of modelling in the optimisation of trial design. Variation of the milk pH and across postpartum time has also been reported and may affect drug distribution (Morriss et al., 1986; Ansell et al., 1977; Matheson et al., 1990), however the intra-individual and inter-individual variability likely outweighs the variability caused by the maturation of breastmilk production (see Supplementary Figure S4).

Breastmilk pH data were available in 15 publications, however this was presented as paired in only half of these (See Supplementary Table S2). Full interpretation of these data is challenging however, when the other parameters which are predicted to affect the M/P ratio are not simultaneously presented. Figure 6 showed that for a number of drugs, the predicted range of M/P ratios was >30-fold, however this is much larger than the observed range of M/P ratios. It is unclear whether this mismatch is due to our incomplete understanding of the effect of pH on M/P ratios (as discussed above), or because the range of milk pH’s collected in each study was more narrow than the physiological range simulated. It is acknowledged that studies quantifying all these milk parameters present a considerable technical challenge, particularly given that milk pH changes with storage time and this likely has knock-on effects for protein binding and lipase activity, which subsequently affects Crt and lipid partitioning (Erickson et al., 2013; Van Den Berg, 1961; Vázquez-Román et al., 2018). It should, however, be feasible to measure milk pH and Crt in milk aliquots prior to storage in order to pair these data with drug concentrations at a later timepoint.

In addition to the quantity of robust data available, another limiting factor in model development may be the quality of data available. To increase the reliability of observed M/P ratios only medications with milk AUC were included, however concerns do still exist. Many of the studies included within the datasets used were conducted over 30 years ago, at a time when analytical methods were less reliable than they are now. The development of deuterated internal standards has significantly improved our ability to quantify drug concentrations within the breastmilk, and our understanding and awareness of drug retention in labware has increased greatly (Alshogran et al., 2024; Lopes et al., 2016). With that in mind, it is noteworthy that the Bartels-LogD model predicted 17/18 non-transporter-mediated drugs within 2-fold of the observed data (across a pH range of 6.63–7.62) when considering publications from 1990 onwards (100% within 3-fold). In contrast, 10/24 drugs published prior to 1990 were mis-predicted, with 9 of these being over predictions. It is plausible to suggest that overpredictions may be more likely to occur due to the reduced reliability of drug extraction and the higher likelihood of sample degradation in older studies. As more clinical data collected in the modern era become available, it is likely that more reliable models may be developed.

Koshimichi et al. developed a semi-mechanistic model for the prediction of M/P ratios using linear regressions to predict the effect of physicochemical properties on secretion (CL_sec) and reuptake (CL_re) into the breastmilk (Koshimichi et al., 2011), with this model being incorporated into the PKsim software to predict infant exposure to 10 medications (80% were predicted within 2-fold; (Nauwelaerts et al., 2023)). Separate prediction of CL_sec and CL_re results in time-dependent differences in the predicted M/P ratio, thus offering an advantage over the models included in the present analysis which calculate a steady state M/P ratio only. The semi-mechanistic model predicted the M/P ratio of 72% of drugs within 3-fold of the observed data, compared with 71%–84% for the Bartels-LogD model (from Figure 6). More recently, an IVIVE model has been developed which incorporates an optimised efflux ratio to improve prediction of transport-mediated drugs (Yang et al., 2022). Across a dataset of 162 drugs (48 of were mediated by passive diffusion only), this model was shown to outperform other models with 66% of 162 drugs being predicted within 2-fold of the observed data compared with 51%, 44% and 43% for the Phase Distribution, Log-Transformed and Koshimichi models, respectively (Yang et al., 2022). This compares with 51%–72% for the Bartels-LogD model.

Overall, these data suggest that the Bartels-LogD model presented here may offer an improvement on existing mathematical models for predicting steady state M/P ratios and performs comparably to other, more mechanistic models published recently. The range of predicted milk-to-plasma ratios presented in this study for each drug demonstrates that consideration of the biological heterogeneity of breastmilk is important for model development and validation. The paucity of data supporting the physiological mechanisms underlying lipid partitioning, protein binding and the effect of ionisation therefore limits the development of models. Although there are technical challenges associated with clinical lactation studies which have hindered this, it is expected that increasing advances in bioanalysis and continuous effort to improve maternal and infant safety will lead to further improvements in model development.