Team:Lund/Model

iGEM Lund 2019

Mathematical Models

We at iGEM Lund have developed mathematical models to see if toxic metal exposure can be treated with a daily intake of our gut bacteria, PRODEACC. This animation shows the of lead accumulation over 8 years when consuming 11 µg lead per day, with and without the use of PRODEACC.

Predicting the Long Term Health Benefits of using PRODEACC

To better understand the long term health benefits of PRODEACC consumption, mathematical models of transformed E. coli Nissle 1917 colonies were developed. Unknown parameters were estimated using parameter optimization against data collected from in-vitro experiments. Once the models could predict the toxic metal bioaccumulation in a single colony, they could be integrated into existing models describing how lead and arsenic are bio-accumulated in the human body. By comparing the output of the modified body accumulation models to the output of the non-modified model, the long term effects of PRODEACC could be assessed. The modeling team has worked in conjunction with the lab group to reach this goal. In the beginning, it was planned to study the deaccumulating effects of PRODEACC on worms in order to assess the long term benefits. This idéa was discarded in favor of using computer models. The computer model has, therefore, had a huge impact on the lab work by alleviating this opportunity cost.


To better understand this text it’s recommended to read the project description before continuing.

Bioaccumulation in a Bacteria Colony

To model the growth, protein production and ion accumulation of a single colony, a system of seven differential equations was devised. Equation (1) and (2) describes the colony as a whole and the concentration of ions in its growth medium, whereas equation (3) - (7) describes the changes inside the average individual cell. A graphical representation of the equation can be found in animation 1. Short descriptions of the time-dependent variables can be found in table 1. Descriptions and approximate values of the parameters can be found in table 2. The initial parameters for the initial conditions can be found in table 3.

Overview of the Bacteria Colony Model

Animation 1: A graphical representation of the rate of change of the ions and proteins described by equation (1) and (3)-(7). (*For the accutal rate, see equation 1.)

Derivation of Equations

The model was devised to simulate the colony in a number of conditions, with only the values of the parameters needed to be changed. For example, for the growth rate kgrowth would most likely be lower in a gut environment compared to an in-vitro environment with good pH and no competition from other microbes.


Equation (1) describes the rate of change of the total number of living bacteria cells in the colony. The cells are assumed to divide at a constant rate kgrowth after the lag phase is over. However, the exponential growth will not continue forever. At some point, the lack of energy will force the colony to stop dividing and will reach a stationary phase. Because of this, the term −kgrowth&times N 2/Nmax was introduced making the growth rate stagnate once N reaches Nmax [1].


Equation (2) describes the rate of change of the total number of lead or arsenic ions in the gut (or growth medium). The rate of accumulation by the colony should be a function of the number of living cells N, as well as the number of transport proteins per cell. The working mechanisms for the lead/iron transporter pbrT are unknown [2]. The Michaelis–Menten equation was used in to represent the transport rate with the rational that a higher concentration of ions outside the cells would increase the probability of an ion finding its way through a transporter, but that the increment would give diminishing returns one more and more transporters work at their highest capacity [3].


Equation (3) describes the change in the number of mRNAs coding for the accumulation proteins and the transport proteins. The transcription rate tsc is assumed to be constant once the promoter has been induced. Then the mRNA will degrade att the kinetic rate kRNA-deg.


Equation (4) describes the rate of change in the number of transport proteins per cell. The proteins are translated from the mRNAs at a kinetic rate ktsl. The proteins are assumed to degrade at the kinetic rate kpro-deg. The kinetic rate -dN/dT/N was introduced to take into account the effective degradation due to cell growth. Every time a cell divides the average number of proteins in one cell is halved. These kinetics apply the other proteins as well.


Equation (5) describes the rate change in the number of free accumulation proteins. Whenever an accumulation protein binds to an ion, the accumulation protein is no longer free. This rate is proportional to the number of accumulation proteins APfree and the number of free ions in cytoplasm C. The kinetic rate kcapture was introduced to scale the rate of capture.


Equation (6) describes the rate of change in the number of occupied accumulation proteins, the proteins bound to an ion. The kinetics are similar to that of equation (5).


Equation (7) describes the rate of change in the number of free ions in the cytoplasm of one bacteria. Ions go from the gut into the cytoplasm, so the equation has the term from equation (2) but with different signs. Once the ions are bound to accumulation proteins, they are no longer free and are removed from the compartment with the same rate as it is removed from the free accumulation proteins compartment from equation (5). Whenever an occupied accumulation protein is degraded, the ion bound to it is released back to the cytoplasm, hence the positive term proportional to the number of occupied accumulation proteins. Last but not least, a negative term to take into account the fact that the number of ions in cytoplasm gets halved every time the bacteria divides.


\frac{dN}{dt} = k_{growth}\times N − k_{growth}\times\frac{N^2}{N_{max}} \frac{dG}{dt} = −\frac{V_{max}}{1+k_m/ \frac{G}{v_{gut}}}\times TP \times N \frac{d(mRNA)}{dt} = tsc − k_{RNA−deg} \times mRNA \frac{d(TP)}{dt} = k_{tsl} \times mRNA − \frac{dN}{dt} \times \frac{TP}{N} − k_{pro−deg} \times TP \frac{d(AP_{free})}{dt} = k_{tsl} \times mRNA − \frac{dN}{dt} \times \frac{AP_{free}}{N} − k_{pro−deg} \times AP_{free} − k_{capture} \times C \times AP_{free} \frac{d(AP_{occupied})}{dt} = k_{capture} \times AP_{free} \times C − \frac{dN}{dt} \times \frac{AP_{occupied}}{N} − k_{pro−deg} \times AP_{occupied} \frac{dC}{dt}= \frac{V_{max}}{1+k_m/ \frac{G}{v_{gut}}} \times TP − k_{capture} \times C \times AP_{free} + k_{pro−deg} \times AP_{occupied} − \frac{dN}{dt} \times \frac{C}{N}

Varible Description
N The total number of living cells in the growth medium.
G The total number of ions in the growth medium.
mRNA The number of mRNAs coding for the ion-transport protein and the ions accumulation protein.
TP The number of transport proteins expressed in a single cell.
APfree The number of free accumulation protein expressed inside a single cell.
APoccupied The number of accumulation protein expressed a single cell, that are bound to an ion.
C The number of free ions inside the cytoplasm of a single cell.


Parameter Description Estimated Value Unit Notes and refereces
Vcell Volume of E.coli cell. 1.5e-10 Litre Allen C. S. Yu et al 2014
Vgut Volume of gut 1.2 Litre Eric Dede 2018
tsc Transcripton rate 900 mRNA/h Genes VIII, Lewin 2004. (Original unit is mRNA/min)
ktsl Translation rate /h Ehrenberg and Kurland, 1988
kRNA-deg RNA degradation rate 8.3177 /h Bionumbers.org, calculated from mRNA’s half life of e.coli which is about 5 minutes.
kpro-deg Protein degradation rate 1.0397 /h Bionumbers.org, calculated from mean protein’s half life in yeast which is about 40 minutes.
kgrowth (in vitro) Cell growth in vitro 1.3862 /h Bionumbers.org Cell doubles every half an hour in vitro
kgrowth (in gut) Cell growth in the gut 0.0173 /h Bionumbers.org Cell doubles every 40 hours in the gut.
Nmax Maximum number of bacteria 10e13 - Thursby E. 2017 Total number of bacteria in human GI tract exceeds 1e14, estimate that around 10% can be E. coli Nissle.
Vmax Maximum transport rate per transport protein 0.35676 ions/TP/h R J Bergeron and W R Weimar The parameter is for another bacteria and another transport protein. Original unit is pg-ion/min/transport protein.
km Ion concentration in the medium, where the rate of transport is half of the maximum. 2.3486e14 ions/m3 R J Bergeron and W R Weimar The parameter is for another bacteria and another transport protein. Original unit is microM.
kcapture Accumulation protein capture rate 0.1 /h Estimated rate since this parameter doesn’t affect the current model


Parameter Description Comments
G0 Number of ions in the gut at start. Should be known by the user when using the model.
N0 Number of cells in the gut/medium at start. Can be set manually or derived from parameter optimization.
tinduction Time of promoter induction Can be set manually.
tlag Lag-phase duration Can be set manually or derived from parameter optimization.

Parameter Optimization

Values for the parameters that couldn't be found in the literature were estimated using data acquired from experiments. The model could be fitted to the data with the help of optimization algorithms. The acquired parameter values were the ones that yielded the lowest cost from a linear loss function. To find the lowest cost adaptive particle swarm algorithm was used[4]. Particle swarm optimization has been used with good efficiency for other non-linear ODE-parameter estimation for biological systems, where there's a big solution space [5].


Particle swarm optimization (PSO) is an optimization algorithm inspired by the movement of a group of organisms, for example, a flock of birds or a shoal of fish. An upper and a lower bound to the search space is given and the algorithm initializes by uniformly distributing the particles/organisms in the space. The solution is calculated for each particle. Each particle then moves in the search space in the direction of the vector sum of is own previous best solution and the group’s best solution. The length at which the particles move varies from different flavors of the algorithm [5].


Particle Swarm Optimization

Animation 2: Animation showing how particles navigate the solution space and converging on the global minimum. Courtesy: Ephramac, CC BY-SA 4.0

Data acquired from wetlab was used in order to derive the correct values for the colony model parameters. It’s important to keep in mind that the parameters are different in different conditions. Which parameter values that could be acquired from which measurement can be seen in table 4.


Experiment Parameters
CFU and OD kgrowth, Nmax, N0, and tlag
SDS-page ktsl
Metal ion concentration in growth medium. ICP-MS was used in this case. Vmax and km

Sensitivity Analysis

When having a complex system of inputs and outputs, it’s difficult to know which inputs have more effects on the outputs. Insight on how to adjust the system to a preferred output can be acquired via sensitivity analysis, in this case how to improve the probiotic’s ability to accumulate metals. Several parameters in the colony model have the potential to be critical, such as kcapture (the capture rate of accumulation protein, which does not have an importance in the current state of this model), km and Vmax (the rate constants for transport protein), or kgrowth and Nmax (the growth rate and the maximum concentration of the bacteria colony), etc. So knowing the importance of these parameters for the de-accumulating effect is also knowing which part of the bacteria colony should be focused on for future improvement. It can either the accumulation protein, the transport protein or the survival of bacteria itself. In order to do that, sensitivity analysis is carried out on the model.


In particular, a global sensitivity analysis (GSA) method called Sobol is used, which is a variance-based analysis. GSA allows us to vary all the parameters simultaneously to evaluate their relative contribution to the model outputs as well as their interaction with each other [6]. Unlike local sensitivity analysis, which only allows for the evaluation of one parameter at a time. The steps for Sobol analysis are described in the figure below.


Sobol Sensitivity Analysis

Figure 1: Pre-Sobol and Sobol sensitivity analysis steps [6]

Before doing sensitivity analysis on parameters, it’s critical to estimate the sensible ranges of them since the analysis would inform the importance of parameters within the ranges. Range estimation is described in table 5. Furthermore, 5000 samples are taken for each parameter within its range. The higher the number of samples is, the more accurate the analysis is. 5000 seems to be a good number of samples for each for a complex system [6].

Parameter Range Range
Nmax 1e10-1e14 The lower bound is the intended initial population of one bacteria colony. The upper bound is based on the fact that there are around 1014 bacteria living in the human gut.
tsc 120-1000 Both lower and upper bound are estimated from the available literature[16].
Vmax 1e-2-1e2 The lower bound is based on literature value used in the parameter table above. Since it is uncertain what the upper bound is, an estimation of 4 order of magnitude bigger than the lower bound is used.
km 1e-5-1e14 The upper bound is based on literature value used in the parameter table above. The lower bound is estimated to be a small number bigger than 0 since it is known that km can’t be negative.
kRNA-deg, ktsl, kpro-deg x/10-x*10 Where x is the value used in the parameter table above. These values are estimated from literature, ranging them one order of magnitude lower and higher seems to be sensible given their values.

Bioaccumulation in the human body

The human body uptake models show how much toxic metal accumulates in different organs of the body. Each of them consists of a system of differential equations, where each differential equation describes the rate of change of the amount of toxic metal in one compartment. A compartment can be either an organ or an excretion pool. We implemented the body uptake models, for arsenic and lead, based on a paper by Eric Dede, published in 2018 [7].

Arsenic

The arsenic body uptake model is adopted from El-Masri and Kenyon (2008) [8] with some modification. In the model, there are 12 compartments that are taken into consideration as seen in figure 1. The figure describes the dynamics between the compartments with first-order rate constants (The constants are left out in this figure for clarity but can be found in Dede 2018. A first-order reaction has a rate proportional to the concentration of one reactant [9], meaning the higher the concentration of a toxic metal is, the faster it is being transferred to other compartments.


Human Body Uptake Model for Arsenic

Figure 2: Overview of compartements and their relations in the arsenic body uptake model.

The model also takes into account that arsenic exists in different forms in the human body. They can be inorganic ions, either As(III) or As(V), or organic forms like monomethyl arsenic acid (MMA) and dimethyl arsenic acid (DMA). There are reactions taking place in the body such that one form of arsenic is converted to another so extra terms are added to the differential equations for that. Since the model contains 12 compartments and 4 forms of arsenic, there’s a total of 48 equations in the system of differential equations, 12 for each of the form of arsenic.


The full system of differential equations including the rate constants can be found in the supplementary materials of Dede et al [7]. To make sure there were no mistakes in the implementation, results were compared to examples in the original paper.

Lead

The body uptake model for lead was originally developed during the 1990s[10][11][12], it was later modified in 2018[7]. There are ten compartments, seven of which are organs while the other three are excretion pools, as seen in figure 3.


Human Body Uptake Model for Lead

Figure 3: Overview of compartements and their relations in non-modfied the lead body uptake model.

To make sure there were no mistakes in the implementation, results were compared to examples in the original paper, when given the same input.


To illustrate how lead accumulates in different organs over many years, an animation was created using data given by the model, see animation 2. This animation illustrates the accumulation when the daily intake is 11.4 µg, over 8 years. After 8 years the accumulation reaches a steady-state for this dose.


Lead Accumulation in the Human Body

Animation 3: The lead accumulation according to the modified body uptake model. The animation shows bioaccumulation over 8 years, when the daily dose was 11.4 µg.

Integration of Colony Model into Human Body Model

The colony models could be integrated into the human body uptake models by adding the colony as a compartment connected to the gut compartment. This way the bacteria colony will compete against the liver when absorbing the metals from the gut. Some assumptions are made for the integration:

  1. The bacteria colony in gut gets excreted out daily.
  2. New bacteria colony is consumed daily.

Arsenic

The colony model was integrated by adding the colony model as a compartment resulting in a model that can be overview in figure 4.


Integrated Human Body Uptake Model for Arsenic

Figure 4: An overview of the bacteria colony integrated human body uptake model for arseic.

Lead

To integrate the colony uptake model into the body uptake it first had to be modified. The original-, modified- and colony integrated versions of the GI-tract-liver dynamics can be seen in figure 5.

Figure 5: Overview of the original liver-gut compartment, the modified liver-gut compartement, and the colony model integrated liver-gut compartement.(*For the actual rate, see equation 1.)


In the original version, the oral lead intake was described as a constant rate IRgi. The rate being absorbed by the liver was described as a fraction of the daily intake. Which means it tells how much of lead is being absorbed daily but not how fast. To integrate the colony uptake model, these dynamics had to be replaced by first-order kinetic rates, similar to the rates found in the body uptake model for arsenic. The daily intake IRgi was then set as an initial condition for the GI tract differential equation. The model could then be solved piecewise, each time with the initial condition G(0) = IRgi. The original (8) and the modified (9) equation for the gut-compartement can be seen below.


\frac{dG}{dt} = IR_{gi} - A_{gi} \times IR_{gi} - (1-A_{gi}) \times IR_{gi} = 0 \frac{dG}{dt} = -k_{GF} \times IR_{gi} - k_{GB} \times G

The values for kGF and kGL could not be found in the literature. To find the values for these rates, while preserving the model's predictability, the parameter optimization, as described earlier, was used.


After the modification, the colony uptake model could be integrated into the lead uptake model, just as in the case with arsenic, resulting in a model that can be overviewed in figure 6.


Integrated Human Body Uptake Model for Lead

Figure 6: An overview of the bacteria colony integrated human body uptake model for lead.

Results

In this section the results from the parameter optimization, the sensitivity analysis, and the model integration are presented.

Parameter Values

The results from the parameter optimization for the colony model can be seen in figure 7, 8 and 9. The values of the parameters can be found in table 5.


Parameter Value Comment
kgrowth 2.31 This is the in-vitro growth rate. It was not used in the integrated model since the growth in the GI-tract is much lower. This value was used during the parameter optimization against levels of concentration of lead.
Nmax 6.12e8 This was used in the parameter optimization but not in the integrated model.
tlag 3.00 This was used in the parameter optimization but not in the integrated model.
Vmax (Lead) 40.9 This was used in the integrated model.
km (Lead) 1.00e4 This was used in the integrated model.
Vmax (Arsenic) 60.0 This was used in the integrated model.
km (Arsenic) 2.60e-2 This was used in the integrated model.
kGB (Lead Body Uptake) 0.44 This was used in the integrated model.
kGF (Lead Body Uptake) 7.02 This was used in the integrated model.

Sensitivity Analysis

The total order for the sobol sensitivity analysis in relation to the outside ion concentration after 3 hours can be seen in figure 10.


Sobol Sensitivity Total Order with respect to Accumulation

Figure 10: Total order sobol sensitivity in relation to the outside ion concentration after 3 hours.

Accumulation in the Human Body

The amount of accumulated arsenic in various organs, with and without the probiotic, can be seen in figure 11 and 12. Figure 11 shows the accumulation after 30 days with a daily intake of 15 µg, whle figure 12 shows the accumulation after 30 days with a daily intake of 50 µg, which is the amount that the general population takes in [15]. The average relative reduction for these scenarios was 83% and 25% respectively.


Accumulation of Arsenic after 30 days when the Daily Dose is 15 µg.

Figure 11: The modeled accumulated arsenic in the kidneys, the liver, the lungs and in the skin, after 30 days with a daily intake of 15 µg per day, with and without the probiotic.


Accumulation of Arsenic after 30 days when the Daily Dose is 50 µg.

Figure 12: The modeled accumulated arsenic in the kidneys, the liver, the lungs and in the skin, after 30 days with a daily intake of 50 µg per day, with and without the probiotic.


Figure 13 shows the number of accumulated lead ions, according to the model, in various organs after 8 years when given a daily intake of 114 µg per day of lead. On average, the relative decreasement when using the integrated models was 8 % for this daily intake.

Accumulation of Lead after 8 years when the Daily Dose is 114 µg.

Figure 13: The modeled accumulated lead in the blood, the bones, the kidnes, the liver, in the poorly- and well-perfused tissue, after 8 years with a daily intake of 114 µg per day, with and without the probiotic.


Figure 14 shows the number of accumulated lead ions, according to the model, in various organs after 8 years when given a daily intake of 11.4 µg per day of lead. On average, the relative decreasement when using the integrated models was 81 %. An animiation illustrating this accumulation was created using the output of the integrated model. See animation 4.


Accumulation of Lead after 8 years when the Daily Dose is 11.4 µg.

Figure 14: The modeled accumulated lead in the blood, the bones, the kidnes, the liver, in the poorly- and well-perfused tissue, after 8 years with a daily intake of 11.4 µg per day, with and without the probiotic.

Animation 4: The accumulation of lead when using the probiotic, when the daily dose is 11.4úg. See animation 3 for comparison.

Discussion

In this section, we will discuss the results presented in the previous section.

Fitted Parameters and Sensitivity Analysis

It is difficult to assess the validity of the parameters in the colon model. More data would need to be collected to know for sure that the parameters are valid. We think it is highly unlikely that km, for arsenic and lead, is correct since it was found by only fitting the model to a single concentration, and since km is the concentration at which the transport rate is half of Vmax. To assess km, one would need many different concentrations when fitting. This would be possible with more, and better, data from the wet lab, but since we outsource the MS-ICP testing to Eurofins, we could not afford more samples. Howecer, from the sensitivity analysis (figure 11) it's evident that the km has a very low effect on the outcome. Therefore, the low uncertainty of the exact value should not have a significant effect on the model's predictability.


The values for Vmax was found to be about two orders of magnitude over the one found for a different transport protein in Ehrenberg and Kurland (1988). We have not been able to judge if our values are beyond reason or not


We have strong confidence in the parameters fitted for the human body uptake model for lead. We tried to reach the authors of the original model, and other experts, to comment on our modification, but without success.


Reduced Accumulation in Human Body

From results, it is evident that the models suggest a significant relative decrease of accumulated lead when daily intake is around 11.4 µg of lead (figure 11), or 15 µg of arsenic (figure 14).


A daily intake of 11.4 µg is far above what is recommended for children by the American Food and Drug Administration (3 µg)[13], but far lower than 114 µg per day, which was observed to be the average daily intake in Kanpur, India in 2005 [14]. For arsenic ic the daily average intake in the US is 50 µg. Since the model did not suggest a big accumulation reduction at these concentrations, we don't believe that this particular probiotic would be useful as a remedy against arsenic exposure. We also do not believe that it would work well against lead at high lead concentrations, for the same reasons (figure 13).


We think that it could be argued that there's a moderate possibility that a probiotic like PRODEACC could be usefull for humans that are moderatly exposed to lead, because the result in figure 14.

Future Improvements

The model we’ve developed can be improved in a number of ways. Here are some that we’ve acknowledged:


The colony model doesn’t take into account bacteria death due to the high concentration of metal in the cytoplasm. This would require another parameter describing the thresh-hold concentration of metal in cytoplasm for causing cell death.


With more data about metal concentration decrease, one could estimate the accumulation parameters with better precision. It would also make it easier to say how good the model represents the features of the colony. It is likely that the value for km is not accurate since the parameter optimization only was done for a single concentration of ions. With more samples from a variety of concentrations, a more accurate value for km would be yielded.


With data from the SDS-page one could estimate the number of expressed protein at different time points and then fit the model to that data to better represent the reality of the protein expression in this colony.


The colony uptake model is developed to from with the lead plasmid in mind, and applied to both lead and arsenic, despite the lack of overexpressed transport proteins in the arsenic plasmid. It, therefore, gives a better representation of the lead plasmid compared to the arsenic plasmid. With more time one could modify the model to better fit the arsenic vector.


The models are implemented in Julia 1.2, and all codes are available on Github for anyone to use or build upon.

Conclusion

Judging from the results, we believe that probiotics, genetically modified or otherwise, are an underappreciated tool for helping the body to get rid of toxic matter. More work needs to be done to evaluate if these particular GM-bacterias can be used in the future as toxic metal remediation. We strongly believe in the methodology used, although the execution can be much improved.

References:

  • [1] Stewart, James; Clegg, Daniel (2012). Brief Applied Calculus. (p. 358)
  • [2] Jaroslawiecka, A., & Piotrowska-Seget, Z. (2013). Lead resistance in micro-organisms. Microbiology, 160(Pt_1), 12–25. doi: 10.1099/mic.0.070284-0
  • [3] Lodish, H. (2002). Molecular cell biology. New York, NY: Freeman. (p. 249)
  • [4] Zhan, Zhang, and Chung. Adaptive particle swarm optimization, IEEE Transactions on Systems, Man, and Cybernetics, Part B: CyberneticsVolume 39, Issue 6, 2009, Pages 1362-1381 (2009)
  • [5] Akman, D., Akman, O., & Schaefer, E. (2018). Parameter Estimation in Ordinary Differential Equations Modeling via Particle Swarm Optimization. Journal of Applied Mathematics, 2018, 1–9. doi: 10.1155/2018/9160793
  • [6] Zhang, X. Y., Trame, M. N., Lesko, L. J., & Schmidt, S. (2015). Sobol Sensitivity Analysis: A Tool to Guide the Development and Evaluation of Systems Pharmacology Models. CPT: pharmacometrics & systems pharmacology, 4(2), 69–79. doi:10.1002/psp4.
  • 6
  • [7] Dede, Eric, et al. "Physiologically-based pharmacokinetic and toxicokinetic models for estimating human exposure to five toxic elements through oral ingestion." Environmental toxicology and pharmacology 57 (2018): 104-114
  • [8] El-Masri, H. A and Kenyon, E. M (2008). Development of a human physiologically based pharmacokinetic (PBPK) model for inorganic arsenic and its mono- and di-methylated metabolites. Journal of Pharmacokinetics and Pharmacodynamics, 35(1), 31–68. doi: 10.1007/s10928-007-9075-z
  • [9] Libretexts. (2019, June 5). The Rate Law. Retrieved from https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Supplemental_Modules_(Physical_and_Theoretical_Chemistry)/Kinetics/Rate_Laws/The_Rate_Law
  • [10] O’Flaherty, E. J. (1991). Physiologically based models for bone-seeking elements. Toxicology and Applied Pharmacology, 111(2), 313–331. doi: 10.1016/0041-008x(91)90033-b
  • [11] O'Flaherty E. J. (1998). A physiologically based kinetic model for lead in children and adults. Environmental health perspectives, 106 Suppl 6(Suppl 6), 1495–1503. doi:10.1289/ehp.98106s61495
  • [12] Fleming, D. E., Chettle, D. R., Webber, C. E., and O’Flaherty, E. J. (1999). The OFlaherty Model of Lead Kinetics: An Evaluation Using Data from a Lead Smelter Population. Toxicology and Applied Pharmacology, 161(1), 100–109. doi: 10.1006/taap.1999.8790
  • [13] Center for Food Safety and Applied Nutrition. (n.d.). Lead in Food, Foodwares, and Dietary Supplements. Retrieved October 19, 2019, from https://www.fda.gov/food/metals/lead-food-foodwares-and-dietary-supplements.
  • [14] Sharma, M., Maheshwari, M., & Morisawa, S. (2005). Dietary and Inhalation Intake of Lead and Estimation of Blood Lead Levels in Adults and Children in Kanpur, India. Risk Analysis, 25(6), 1573–1588. doi: 10.1111/j.1539-6924.2005.00683.x
  • [15] Agency for Toxic Substances and Disease. (2007, August 1). Toxic Substances Portal - Arsenic. Retrieved October 19, 2019, from https://www.atsdr.cdc.gov/phs/phs.asp?id=18&tid=3.
  • [16] https://bionumbers.hms.harvard.edu/bionumber.aspx?id=111997&ver=5