Structure of the PBTK modelCore model structure

The core model structure of our Al PBTK model is identical to our previous model from Hethey et al. (2021). It is comprised of ordinary differential equations (ODEs) describing the rate of change of Al in blood and tissues (bone, brain, liver, kidney, spleen, muscle, rest of body), combined with dosing compartments and an elimination compartment (urine), with a physiological parametrization typical for PBTK models (Jones and Rowland-Yeo 2013). In blood, two different Al species are modelled, namely “Mix” (equilibrium of Al salts) and “addCit” (excess Al citrate), with a first-order transfer term describing equilibration of Al salts after i.v. dosing of Al citrate. Distribution into tissue is modelled without distinguishing between Al species, due to a lack of sufficient quantitative understanding of the uptake mechanisms. \(\textrm\) redistribution from tissue to blood is divided into addCit and Mix according to the current addCit:Mix ratio. For all tissues except bone, the uptake and release rate constants are parametrised as follows:

$$\begin k_\mathrm }2}} = \frac\cdot Q_\textrm}}, \qquad k_\mathrm }2}} = \frac}\cdot V_\textrm}, \end$$

(1)

with tissue blood flow \(Q_\textrm\), total blood volume \(V_\textrm\) and tissue volume \(V_\textrm\) being physiological parameters and the uptake and retention coefficients \(I_\textrm\) and \(K_\textrm\) estimated parameters, with \(I_\textrm\) empirically modelling tissue uptake processes not fully understood mechanistically.

The data basis for parameter estimation was a comprehensive collection of \(^\textrm\) tracer studies from the literature (Weisser et al. 2017; Hethey et al. 2021), focussing on rats and humans and Al chloride and Al citrate as the most frequently studied species and (soluble) Al salts, respectively. The choice of tissues was also based on data availability. Physiological and estimated parameters were assumed to be tissue-specific, except that visceral organs (liver, spleen, kidney) were assumed to share uptake and retention coefficients (due to common physiological properties and data sparsity). More details on the above processes and assumptions can be found in Hethey et al. (2021).

The newly implemented modules for distribution into bone and renal clearance, the implementation of various dosing routes and physiological growth are discussed in separate sections below. The extended model structure is illustrated in Fig. 1, also indicating updates with respect to Hethey et al. (2021). The full system of ODEs is provided in the Supplementary Material, Sect. S1.

Fig. 1

Structure of the Al PBTK model. Each compartment represents an ordinary differential equation, straight arrows indicate an exchange of \(\textrm\) between compartments and jagged arrows are dosing sites. All exchange processes were modelled as first-order, except that absorption of Al adjuvants was modelled as a zero-order process. Symbols written on top of arrows are model parameters, with (non-estimated) physiological parameters in parentheses. Modules/parametrizations which have changed or are newly added compared to Hethey et al. (2021) are highlighted in red. Abbreviations: AlCit, aluminium citrate; AlChl, aluminium chloride; PN, parenteral nutrition; i.v., intravenous; p.o., per oral; s.c., subcutaneous; i.m. intramuscular

Dynamic bone module based on calcium kinetics

In our previous model (Hethey et al. 2021), bone uptake and release rates were parametrised using Eq. 1, like any other organ. While this parametrisation reflected the effects of net bone growth, it did not capture the known changes in bone turnover during childhood and adolescence.

Similarly to other models for bone-seeking elements like lead (O’Flaherty 1993), and strontium (Pertinez et al. 2013), we therefore linked \(\textrm\) uptake into and release from bone to that of Ca, thereby leveraging the detailed knowledge on changes in Ca kinetics during growth. We recently developed a consolidated age- and sex-dependent description of Ca deposition and release rates (\(v^\text _\) and \(v^\text _\), respectively) into bone in humans, which is moreover consistent with the concurrent changes in body size (weight, bone mass) (Hartung et al. 2024). To be usable in the Al PBTK model, these rates first needed to be transformed into first-order uptake and release constants via

$$\begin k_\text ^\text = \frac_}}\cdot C^\text _\textrm};\qquad k_\text ^\text = \frac_}}\cdot C^\text _\textrm}. \end$$

Ca concentration in human blood \(C^\text _\textrm\) is stable over the lifetime, and was fixed was fixed to 2.5 mmol/L based on literature (Sofronescu 2019), while Ca concentration in bone was determined based on bone calcium content and bone density, i.e. \(C^\text _\textrm= f_\text \cdot D_\text \), from ICRP (2002). In rats, we assumed \(C^\text _\textrm\) = 2.25 mmol/L (Granjon et al. 2016) and derived Ca uptake and release constants \(v^\text _\) and \(v^\text _\), as well as calcium concentration in bone \(C^\text _\textrm\), from a rat bone growth model proposed in O’Flaherty (1991a, 1991b). Details are given in Supplementary Material section S2.

Then, first-order \(\textrm\) uptake and release constants were determined via

$$\begin k_\text ^\text = \alpha \cdot k_\text ^\text ;\qquad k_\text ^\text = \beta \cdot k_\text ^\text , \end$$

with species-independent \(\alpha\) and \(\beta\). Since mechanistically, release of Ca and Al from bone are both assumed to be due to bone resorption by osteoclasts, the typical value of \(\beta\) was fixed to 1, meaning that the typical Al release rate is equal to the reference release rate determined for Ca (however, random inter-individual variability around this value was included). Al uptake from blood to bone is assumed proportional to that of Ca (and not equal) since uptake mechanisms differ, similarly to other TK models for bone-seeking elements (e.g. Leggett (1992), Assumption1). In total, this bone module contained three parameters (fixed+random effects for \(\alpha\) and random effect for \(\beta\)), compared to four in the previous model version (fixed+random effects for bone uptake and retention coefficients).

Renal clearance module

Al is almost exclusively excreted via the kidneys, mainly as Al citrate (Priest 2004; Shirley and Lote 2005). The previous model in Hethey et al. (2021) featured an effective renal clearance (CL), which accounted for the Al citrate fraction in blood. While not being a problem for single doses as considered in Hethey et al. (2021), this formulation proved to be inappropriate for combined dosing schedules (e.g. p.o. and i.v. routes in parallel), where it led to implausible (non-additive) interactions. We therefore replaced the effective CL by separate Cit and Mix clearances,

$$\begin \text _\text = \frac\cdot \text _\text \cdot \text }+\text _\text \cdot \text }, \qquad \text _\text = \frac\cdot \text _\text \cdot \text }+\text _\text \cdot \text }, \end$$

with the same ultrafiltrable fractions \(\text _\text = 1\) and \(\text _\text = 0.1\) as in (Hethey et al. 2021), and GFR derived from species-specific models as explained in section 2.1.5. During model development in Hethey et al. (2021), estimation issues were encountered with such a model, but together with other reformulations, a successful estimation was now possible (see Sect. 3.1).

Dosing module

The previous model in Hethey et al. (2021) featured i.v. bolus dosing of Al citrate and chloride (into compartments “addCit” and “Mix”, respectively), as well as p.o. dosing, with a common absorption rate and bioavailability for the two salts and assuming systemic uptake via compartment “Mix”. We extended the dosing module to include other routes as well:

Al exposure via food was modelled using the same bioavailability and absorption rate as for p.o. dosing, but as a continuous uptake rate rather than discrete (meal) events;

For orally dosed Al-based antacids, we used the reported bioavailability for the specific antacid, together with the estimated absorption rate for soluble Al salts;

Al exposure through PN was treated the same way as i.v. dosing of Al chloride, since PN solutions are expected to mainly contain this Al salt;

s.c. and i.m. dosing of Al adjuvants was modelled using zero-order release kinetics from a depot compartment as Al citrate. To this end, injection site release data after s.c. and i.m. administration of different Al adjuvants (Al hydroxide (AH) or Al phosphate (AP), either prepared in situ or as commercial gel formulations (e.g. Alhydrogel, Adju-Phos)) in various animal species (rat, rabbit, monkey) were collected from literature (Flarend et al. 1997; McDougall et al. 2016; Verdier et al. 2005; Weisser et al. 2019, 2020). Based on these data, zero-order absorption constants were estimated for each adjuvant, based on a linear regression with fixed intercept at 100%. These absorption kinetics were used for simulations of Al-containing adjuvants in both rats and humans.

Age-dependent physiology

In Hethey et al. (2021), all considered human physiological and anatomical parameter values were from adults and assumed as static. For children, however, age-dependent physiological changes during Al exposure need to be accounted for, especially when considering slow processes such as Al release from bone. We therefore rendered all physiological and anatomical parameters age-dependent, from (term) newborns to adults, as explained below. For most parameters, literature values for different age groups (newborn, 1 year, 5 years, 10 years, 15 years and adult) were linearly interpolated: body weight, organ weights and bone density were taken from ICRP (2002), and blood flows were derived using a lean body weight/fat mass-based scaling approach proposed in Huisinga et al. (2012). Of note, the definition of the bone compartment was changed compared to Hethey et al. (2021), now representing cartilage-free bone instead of total bone comprising cartilage (see Sec. S3 for details). For Ca-related parameters and GFR, specific models were used. Ca bone deposition and release rates \(v_^\text \) and \(v_^\text \) were chosen as explained in Sec. 2.1.2. GFR (in mL/min) was calculated based on an allometric scaling/maturation approach described in Rhodin et al. (2009),

$$\begin \text = \text _\text \cdot \left( \frac}_\text }\right) ^ \cdot \frac^}^ + \text _^} \end$$

(2)

where lean body weight \(\text \) is predicted based on body weight and height using a model by Janmahasatian et al. (2005), PMA is post-menstrual age (age from birth plus a pregnancy duration of 0.77 y), and with parameters \(\text _\text = \,\mathrm\), \(\text _\text = \,\textrm\), and \(\text _ = \,\textrm\), the latter being the maturation half-time. The age-dependency of all physiological parameters is illustrated in Fig. S4.

For rats, static physiologies were used, as in Hethey et al. (2021). For each simulation, the rat body weight was chosen to be either 250 g or 475 g, with associated detailed physiologies reported in Hethey et al. (2021). For studies with short duration (\(\le\)5 d), the reference rat with the smallest difference to the mean body weight at study onset was chosen. In case of longer study duration, the number of data points generated at time points \(\le\)5 d and >5d were counted. If the majority of data points were \(\le\)5d, the starting weight of the rats were used for allocation. Otherwise, the heavy reference rat (475 g) was used, since rats are growing rapidly and reach their adult weight within a few weeks (e.g. a 300 g Wistar rat has grown to about 400–450 g within 30 days, Janvier (2024)). Rat GFR was predicted based on body weight (BW) using the same approach as described in Hethey et al. (2021), namely a purely allometric relationship

$$\begin \text = \text _\text \cdot \left( \frac}_\text }\right) ^ \end$$

(3)

with \(\text _\text = \,\mathrm\) reported by Davies and Morris (1993) for a rat with body weight \(\text _\text = \,\textrm\).

Simulation details

All numerical simulation results were obtained using R version 4.2.2 (R Core Team 2022), using package mlxR version 4.2.0 (Lavielle 2021) for simulation of the Al PBTK model. For all exposure scenarios, \(n=500\) individuals were simulated and simulation results displayed as median and 90% confidence intervals. Simulations compared to individual data were always matched by age and sex. In contrast, for simulations compared to summary statistics, a mixed male/female population was simulated and the average reported age used. R Markdown scripts to produced all results are available on Zenodo (https://doi.org/10.5281/zenodo.14710873).

Initial conditions

Especially for children below 1 year of age, initial Al levels built up during embryonal development need to be considered. To this end, the literature was screened for reports on normal Al concentrations in healthy newborns. Six references were identified, reporting on one or several organs amongst plasma, bone, brain, liver and kidney (Moreno et al. 1994; Bozynski et al. 1989; Sedman et al. 1985; Litov et al. 1989; Bougle et al. 1992; Hawkins et al. 1994). These reported values were used to identify the parameters of an assumed underlying lognormal distribution of Al concentrations (for details, see Section S6.1). Separately for each organ, lognormal distributions obtained this way from different literature sources were averaged arithmetically; subsequently, another lognormal distribution was fitted to the averaged distribution, yielding the final distribution from which to sample initial values. This process is shown in Fig. S3. For muscle and spleen, no literature values on Al concentrations in newborns could be retrieved. Spleen being a visceral organ like liver and kidney and the Al concentration distributions for the latter two organs being very similar, an average of the liver and kidney distributions profiles was assumed for spleen (consistent with the use of common intrusion and retention coefficients \(I_\textrm,K_\textrm\) in visceral organs in the Al PBTK model, leading to identical steady-state concentrations). For muscle as well as for rest of body, both of which are in quick exchange with plasma, a quasi-steady state with plasma was assumed, based on the parameter estimates reported in Table S2.

We found no data on the correlation structure of initial values in different organs. A plausible assumption is that if the initial value is large in one organ, this is also the case in other organs. We used a simple mathematical representation thereof, assuming the distributions in different organs to be perfectly correlated. Algorithmically, this means that, for each simulation, a single standard normally distributed value was drawn, which was then scaled to match mean and variance of each normal distribution underlying the lognormal initial distributions in different tissues, and exponentiated.

For consistency, these initial value distributions were applied for all simulations of \(^\textrm\)-based studies in humans (children and adults), except for Klein et al. (1982); Stoehr et al. (2006) which start at adult age, neglecting baseline exposure due to the long duration and large exposure of the PN solution.

Reference Al intake via food

For simulation scenarios of \(^\textrm\) exposure, the underlying intake by food (depending on age) needs to be included, except during PN. An average realistic Al intake via food used in the simulations was taken from a review of European studies published by EFSA (2008); see Table 1 for an overview and Supplementary Material section S5 for details.

Table 1 Average Al intake via food based on European studies on food consumption (EFSA 2008), which was used in the simulations

Since the estimated bioavailability (F) of 0.168% in our model (see Table S2), representing an average for Al chloride and Al citrate in solution, was in good agreement with commonly accepted values of 0.1% (food) – 0.3% (drinking water) for Al absorption from diet (Tietz et al. 2019; Yokel and McNamara 2001), this estimate was used for exposure from food as well.

Dry versus wet tissue weight

The content of Al measured in bone and other tissues is usually reported as \(\mu\)g/g dry weight or wet weight, depending on the condition of the sample when it was weighed. In the Al PBTK model, all organs including bone are considered as “wet”. Therefore, all literature data reported as dry weight were transformed to wet weight to allow for a meaningful comparison. In addition, bone samples were considered as marrow-free bone tissue (which was mostly, but not always clearly indicated), whereas the model bone compartment comprises cartilage-free skeleton including bone marrow, which also needed to be accounted for when transforming the data. Conversion factors to derive wet weight (model) concentrations from dry weight data in humans are given in Table 2; their derivation from literature values can be found in Supplementary Material section S4.

Table 2 Dry-to-wet weight conversion factors (C), such that 1 \(\upmu \textrm/\textrm\) dry weight \(\hat~C\) \(\upmu \textrm/\textrm\) wet weightLiterature search and summary of validation datasets

The literature was screened for \(\textrm\) concentration data under diverse conditions to evaluate the model’s capacity to correctly predict \(\textrm\) exposure. In contrast to our previously assembled \(^\textrm\) dataset in Hethey et al. (2021), here we also considered \(^\textrm\) data to allow for a more thorough validation. In total, 22 references were selected for model validation (see Table 3), as explained below.

Normal \(^\textrm\) plasma and tissue concentrations in healthy adults allow for an evaluation of the model’s capacity to correctly predict the impact of long-term \(\textrm\) exposure, presumably mainly via food. Only studies explicitly designed to collect normal values (not control groups in non-representative populations such as hospitalized patients) were considered. Reference published more than 30 years ago were excluded due to potential changes in \(\textrm\) exposure over time. We identified 11 references reporting reference \(^\textrm\) concentrations in different tissues in adults. Besides in plasma and toxicologically relevant organs (bone, brain and liver), we also found reference values in kidney, spleen and muscle. These values were collected from both US/Canadian (Krewski et al. 2007; Tang et al. 1999; Kruger et al. 2014; McLachlan et al. 2018) and European populations (Hellström et al. 2005, 2006; House et al. 2012; Exley and Clarkson 2020; Roider and Drasch 1999; Rahil-Khazen et al. 2002; Hongve et al. 1996). Two studies reported results based on wet tissue weight (Roider and Drasch 1999; Rahil-Khazen et al. 2002), the other nine studies in terms of dry tissue weight. The latter were converted to wet weight as described in Sect. 2.2.3.

PN, and especially total (i.e. exclusive) PN, may result in elevated \(^\textrm\) exposure. To identify suitable validation datasets, we performed a literature search in Pubmed and Google scholar using the search terms “aluminium” AND “parenteral nutrition”. The results were kept if they contained data on both Al measurements in plasma or tissues during/directly after PN and the cumulative PN volume administered up to the sampling timepoint(s) including Al content, allowing calculation of the total Al intak