A pharmacokinetic model of lead absorption and calcium competitive dynamics

Lead is a naturally-occurring element. It has been known to man for a long time, and it is one of the longest established poisons. The current consensus is that no level of lead exposure should be deemed “safe”. New evidence regarding the blood levels at which morbidities occur has prompted the CDC to reduce the screening guideline of 10 μg/dl to 2 μg/dl. Measurable cognitive decline (reduced IQ, academic deficits) have been found to occur at levels below 10 μg/dl, especially in children. Knowledge of lead pharmacology allows us to better understand its absorption and metabolization, mechanisms that produce its medical consequences. Based upon an original and very simplified compartmental model of Rabinowitz (1973) with only three major compartments (blood, bone and soft tissue), extensive biophysical models sprouted over the following two decades. However, none of these models have been specifically designed to use new knowledge of lead molecular dynamics to understand its deleterious effects on the brain. We build and analyze a compartmental model of lead pharmacokinetics, focused specifically on addressing neurotoxicity. We use traditional phase space methods, parameter sensitivity analysis and bifurcation theory to study the transitions in the system’s behavior in response to various physiological parameters. We conclude that modeling the complex interaction of lead and calcium along their dynamic trajectory may successfully explain counter-intuitive effects on systemic function and neural behavior which could not be addressed by existing linear models. Our results encourage further efforts towards using nonlinear phenomenology in conjunction with empirically driven system parameters, to obtain a biophysical model able to provide clinical assessments and predictions.

In the 1980s, leaded gasoline was deemed "environmentally unsafe" and forced out of the market place. Regulatory mechanisms limiting the content of lead and its compounds in paint and gasoline were remarkably successful in reducing the prevalence of highly elevated blood lead levels in both adults and children. Studies have shown that the mean blood lead level of persons aged 1 to 74 years dropped 78% (from 0.62 to 0.14 μmol/l) from 1976 to 1991, attributed primarily to the almost complete removal of lead from gasoline and soldered cans 8 .
These regulatory actions lead in the 1990s to the premature hope that the lead contamination problem had been solved. However, despite the constant effort and increasing success in eliminating sources of exposure, lead remains the most important pediatric environmental health problem, with costs associated with lead-related morbidities estimated in the billions of dollars 7 . While cases of lethal intoxication are currently extremely rare, a more major concern consists of the now documented neurodevelopmental effects resulting from children's continuing exposure to low levels of lead. New evidence regarding the blood lead levels at which morbidities occur have been putting pressure on the CDC to reduce the current screening guideline of 10 μg/dl to 2 μg/dl 9 . Measurable cognitive decline (reduced intelligence quotient and academic deficits) have been found to occur at levels below 10 μg/dl. Increased exposure has also been associated with neuropsychiatric abnormalities, including attention deficit hyperactivity disorder and antisocial behavior 7 . See also 10,11 and 12 . Functional imaging studies are beginning to shed some light onto the neural mechanisms of the neurodevelopmental effects of lead. Knowledge of lead pharmacology allows us to better understand its processes of absorption and metabolization and the mechanisms that produce its medical consequences, to strive for ways to prevent contamination, and to facilitate treatment after exposure 13 . The consensus is that no level of lead exposure appears to be "safe, " with some studies even suggesting that the rate of decline in performance is greater at levels below 10 mg/dl than above 10 mg/dl. However, no plausible mechanisms have been identified 7 . A biophysically informed mathematical modeling approach may provide a framework for phrasing some of these unanswered questions, and a path to understanding the subtle effects that occur at these low levels of exposure.

Basic Lead Biokinetics
No biological requirement for lead has ever been demonstrated, and the human body does not metabolize it into other elements. Lead typically enters the body in a few ways: via ingestion (through the digestive system), through breathing (via airways) and in small quantities through skin. If received via the gastrointestinal pathway, the effectiveness of the absorption depends on the individual's food intake prior to exposure, both quantitative (since food consumption decreases absorption of water-soluble lead) and qualitative (due to interactions with other elements in the diet). It is also well-known that efficiency of gastrointestinal absorption of water-soluble lead is age dependent, and substantially higher in children than in adults.
Statistical reports from studies of soft tissue concentrations of lead, in both humans and animals, have changed dramatically over the years, from the age when occupational exposure levels were high [14][15][16][17] to more recent years, characterized by lower exposure 18 . Throughout the downward trends in soft tissue lead levels, autopsy studies provide a basis for describing the relative soft tissue distribution of lead in adults and children. Most of the lead in soft tissue is in liver 14,16,[19][20][21] . In this study, however, we will be primarily interested in the dynamics and effects of lead on neural tissue, which will therefore be studied and discussed as a separate soft tissue compartment.
Approximately 95% of lead in adult tissues, and approximately 70% in children, resides in mineralized tissues such as bone and teeth 14,15 . This reflects changing turnover rates along an individual's lifetime, with a slower turnover of lead in adult bone than in children [14][15][16][17] . The lead deposit in adult bone can act to replenish lead eliminated from blood by excretion, even long after exposure has ended [22][23][24][25][26] . It can also act as a source of lead transfer to the fetus when maternal bone is resorbed for the production of the fetal skeleton [27][28][29][30] .
The main excretion pathway for lead is via kidney clearance; other paths such as via sweat, saliva, hair and fingernails, are by comparison negligible. Mechanisms by which inorganic lead is excreted in urine have not been fully characterized. Such studies have been hampered by the difficulties associated with measuring ultrafilterable lead in plasma and thereby in measuring the rate of glomerular filtration of lead 31 . Measurement of the renal clearance of ultrafilterable lead in plasma indicates that lead undergoes glomerular filtration and net tubular reabsorption 32,33 . Renal clearance of blood lead increases with increasing blood lead concentrations above 25 μg/dL 34 . The mechanism for this has not been elucidated and could involve a shift in the distribution of lead in blood towards a fraction having a higher glomerular filtration rate (e.g., lower molecular weight complex), a capacity-limited mechanism in the tubular reabsorption of lead, or the effects of lead-induced nephrotoxicity on lead reabsorption.
Over the past half a century, since the age of peak lead exposure in the 1960s, models of led dynamics and pharmacokinetics have evolved substantially. Based upon the original and very simplistic compartmental model of Rabinowitz 35 , built in 1973 with only three major compartments, a lot of more extensive biophysical models sprouted over the following two decades, using carefully documented physiological and biological empirical estimations of the inter-compartment rates and other model parameters. Among these, the Lead Metabolism Model of Rabinowitz et al. 36 distinguishes two different soft tissue compartments (deep and shallow soft tissue, each with different lead dynamics), and the Marcus Model 37-39 considers two different bone compartments (cortical and trabecular). The O'Flaherty Model additionally distinguishes well-perfused and poorly-perfused tissues [40][41][42][43][44] . The IEUBK (Integrated Exposure Uptake Biokinetic) model 45,46 set new steps in differentiating between sources of lead contamination, and was used to investigate age effects (particularly in children). Since the 1990s, the Leggett Model 47 has been used as the state of the art compartmental chart for lead pharmacokinetics, with over 20 compartments, and introducing more subtle, nonlinear flow rates to reflect age effects. A recent CDC study 48 draws a thorough comparison between models, emphasizing their strengths, and comparing the efficiency of their risk assessment.
This large body of modeling work has been very useful in understanding the path of lead through the body, and generate predictions of toxic levels based only on a priori knowledge of intake and biophysical characteristics Age dependence and effects in children. The toxic effects of lead have been overwhelmingly observed in children. Some of the health effects that have been associated with lead exposure in children are similar with those observed in adults at higher exposures. However, children's susceptibility and response may qualitatively and quantitatively differ from those encountered in adults, due to their specific physiology and behavior, which can influence both exposure and processing. The effects depend on developmental age (which influences pharmacokinetics and metabolism 49,50 . Typically, the peak of the vulnerability and disruptive effects occurs during critical periods of structural and functional development. Health effects include anemia 51 , renal alterations, impaired metabolism of vitamin D 52,53 , growth retardation, delayed puberty 54 . Exposure to lead during childhood is also well-known to result in neurobehavioral effects that persist into adulthood and may not be evident until a later stage of development (making difficult a correlation, or a precise quantitative assessment of these effects). Delays or impairment in neurological and neurobehavioral development have been noted even at very low doses, and include encelopathy, lower cognitive performance 11,55 , neuropsychiatric disorders such as attention deficit hyperactivity disorder and antisocial behavior 7 .
To start with, behavioral patterns of children can result in higher rates of ingestion of lead (e.g., from soil and dust, both of which are often important environmental depots for lead [56][57][58] . Children also absorb a larger fraction of ingested lead than do adults; thus, children will experience a higher internal lead dose per unit of body mass than adults at similar exposure concentrations [59][60][61] . It was suggested that the gastrointestinal absorption of lead is greatest in infants and young children 61 . There may also be differences in excretion, since infants have lower glomerular filtration rate an inefficient tubular secretion and resorption capacities 49,62 .
While toxicokinetics of lead in children appears to be similar to that in adults, the action of many xenobiotic metabolizing enzymes seem to depend on developmental stage [63][64][65] . This lead to the children's increased susceptibility to toxic effects of lead or to detoxification, although the exact mechanisms of this sensitivity is not completely understood. Children and adults may differ in their capacity to repair damage from the deleterious effects of lead poisoning. However, it is important to note that children also have a longer remaining lifetime in which to express damage from toxic exposure.
Several models of lead pharmacokinetics in children have been developed [43][44][45]47 . We considered these models when incorporating dependence on age in our own mathematical model.

Modeling Methods
In the present study, we construct and analyze a compartmental model of lead pharmacokinetics. We are primarily interested in understanding the neurotoxic effects of lead, and their modulation across development, as described in the previous section). Neural effects of lead contamination are tightly related to the pharmacokinetics of calcium, in ways which have been clearly established at most points of their dynamic trajectory through the organism.
It was shown that the competitive presence of calcium can affect: (1) lead's intestinal absorption, (2) its kinetics between soft tissues; (3) its storage in bones and its mobilization from osseous to non-osseous tissue; (4) its retention versus excretion rates; (5) the toxic response of the body to lead. Early research in rodents has revealed that a lower Ca diet increased their susceptibility to the toxic effects of lead 66 (including order of magnitude higher lead blood levels, anemia, renal problems). A few mechanisms were proposed 67 . First, it was suggested that higher Ca directly decreases Pb absorption by creating competition for its transport through the gastro-intestinal phospholipidic wall 68,69 . Second, it was also found that decreased Ca acts as an inhibitor for the Pb release from the skeleton, thus increasing the body Pb content 70 . Third, it was speculated that the increase in Pb due to low Ca involves the kidney as the action site, but no mechanism was proposed 71 . More recently, the Ca/Pb competitive interaction has been noted more generally, at other interfaces of the system, including the blood/brain barrier (as further described in our Modeling Methods). Since lead competes with calcium on transporters, lead can also replace calcium, with serious consequences on bone formation, kidney function and, most importantly, neural function (where calcium is indispensable for processes like learning and memory).
We therefore aim to study simultaneously lead and calcium kinetics, and their tight interaction all along their dynamic trajectory through the body. In order to do so, we build a middle ground model (see Fig. 1), which follows the traditional idea of Rabinowitz, rather than the newer models (more elaborate on the compartment number and specificity). We will rather aim to retain as much as possible of the model's simplicity and low-dimensional nature, while trying to incorporate and understand the more subtle, calcium-mediated, nonlinear aspects of the transmission between compartments.
Research has shed increasingly more light over the past few decades on the role of other elements (besides Calcium) to modulating the body's absorption and processing of lead, as well as its susceptibility to lead's toxic effects. Together with macroscopic measures and observations, there is increasing knowledge of molecular competition between lead and these other elements (such as Iron, Phosphates, Zinc) when crossing cellular membranes. Our model may be a first step towards building a more general, unified quantitative theory addressing how effectively feeding and nutritional patterns may be used to minimize lead absorption and toxic effects, or to We build upon the original variables used in the Rabinowitz model (blood, soft tissue and bone), and refine them to distinguish the brain and the kidneys as separate soft tissue compartments, as shown in Fig. 2.
We concentrate on lead intake through the digestive tract (and ignore other sources). After ingestion, environmental lead p and calcium c get absorbed into blood through the gastric endothelial cells, via a common/ competitive mechanism. We use the time variable S Pb and S Ca the designate the blood content of lead and calcium, respectively.
From blood, lead passes into soft tissue. In previous models, this compartment has been further separated, based on the aim of the model, into either shallow and deep soft tissue 36 , or into organs or even organ parts 47 . For the needs of our model, we only consider the brain (cognitive function) and the kidneys (excretion) as separate Compartments included in our model. Environmental lead enters the blood stream (primarily via the digestive tract); from blood it crosses into soft tissue (which we further separate into brain, kidneys and other soft tissue), or it deposits in bones, from where it can cross back into the blood stream, under favorable biochemical circumstances. From kidneys it is excreted via urine (other forms of excretion are considered negligible and ignored in this model). Lead in the brain impairs normal neural function.

Figure 2.
Model of coupled lead and calcium kinetics. Compartments: S = blood; O = bone; B = brain; K = kidney; T = other soft tissue. In each compartments, the subscript Pb refers to lead content, and Ca to calcium content in the corresponding compartment. The parameters p and c represent ingestion rates for lead and calcium, respectively. The arrows show the directions of the two elements leaving each compartment X and entering a new compartment Y, through shared molecular transporters, adding to the content of the new compartment Y. The equations that govern the flow between compartments are further explained in the main text.
compartments, and not focus on other processes (such as liver metabolization, considered in previous models). We will refer to the rest as "other soft tissue, " and use the variables T Pb and T Ca to designate lead and calcium content in the other soft tissue, respectively.
We support the idea of the brain as a stand alone variable for two reasons. First, we expect that lead transport to occur differently between the blood and the brain than between the blood and any other soft tissue, due to the additional protection provided by the blood-brain barrier (BBB, which we incorporate in the model). Secondly, as explained in the introduction, our efforts are invested in understanding the deleterious effects of lead toxicity on brain function. Absent the massive lead exposure associated to the 1960s-70s, recent research has focused on investigating the effects on cognition and behavior of subtle lead exposure (in children in particular). Children have an immature blood-brain barrier 72,73 . Recent studies point out a nonlinear dose-effect dependence, so that relatively small elevations of blood levels of lead in the low level contamination range may produce on cognitive performance effects that are proportionally greater than those observed in conjunction with higher exposure and with more substantially elevated blood levels 10,[74][75][76][77] . No physiological explanation has yet been proposed for this nonlinear effect, which is key in understanding and addressing lead toxicity derived cognitive deficits in children.
We investigate whether a model can shed some light on potential mechanisms that underlie the nonlinearity. We call B Pb and B Ca the lead and calcium levels in brain tissue, the balance of which impact brain function in a variety of ways, among which is the effect of calcium on neurotransmitters (in particular on the function of glutamate).
We also treat the kidneys as a separate module, for two reasons. First, because they are the main gateway for lead out of the body (other means, such as via hair of nail growth, being small enough to be considered negligible). Secondly, it has been noted that lead toxicity may lead to kidney dysfunction, thus producing a broken feedback loop which has the potential to exacerbate both renal disease and toxic effects of lead (by diminishing its excretion). Recent studies are increasingly documenting potential molecular mechanisms of this interaction, making it an ideal candidate for mathematical modeling. A model can be used to further investigate this mechanism, make predictions and optimize medication plans. We call K Pb and K Ca the kidney content of lead and calcium, respectively.
Most of the lead body burden resides in bones, with a slightly larger fraction in adults compared to children 14 . There is insufficient information to determine how similar the lead metabolism in children is to that in adults. The mechanisms of lead deposit to bone, and lead release from bone are clearly tied with the dynamics of calcium, and vary greatly along an individual's life (with kinetics increasing drastically during time windows of high calcium resorption/mobilization, such as pregnancy, or old age). Our model will explore the relationship between lead and calcium dynamics at this level as well. For simplicity, we will do so without distinguishing between different bone layers (as has been done in other existing models [37][38][39][40]43,47 ), and we will consider overall variables O Pb and O Ca for the lead and calcium bone content, respectively.

inter-compartmental Dynamics
It is known that transfer of lead from one compartment to another typically involves saturable molecular mechanisms of lead-calcium competition. Let X Pb and X Ca be the variables describing the content of lead and calcium in a certain compartment X, respectively. The two types of molecules will compete on a fixed number of transporters in a nonlinear fashion depending on the "crowdedness" (i.e., joint number of molecules + X X Pb Ca in the compartment), so that increasing the pool will first increase the flow slowly, than more efficiently, but eventually would saturate due to the competition over the same number of transporters. Then the likelihood of one molecule (of either lead or calcium) to move from X into Y can be expressed as a traditional sigmoidal function: which satisfies the following properties: it is increasing from 0 (as → −∞ x ) towards a saturation value of 1 (as → ∞ x ); it goes through an inflection point and a high sensitivity window whose position and width are modulated by the parameters b and θ. In order to fix our ideas, and illustrate the principles behind the model behavior (with respect to variations in lead/calcium intake, as well as in other physiological parameters of interest), we fixed the sigmoidal parameters to the values = .
b 0 6 and θ = 6 for all compartmental pairs, throughout the analysis. As shown in Appendix A, the behavior of the system is qualitatively robust under variations of θ b ( , ) within a relatively large range, implying that the results and phase transitions we present in the forthcoming sections are not significantly affected by fine tuning of the sensitivity of the molecular transport. The values = .
b 0 6 and θ = 6 were particularly chosen within this range as the values that produced quantitatively the best aligned results with existing empirical information, as further discussed in the Results section.
The probability for a molecule that crosses from X to Y to be a molecule from X Pb can be simply modeled as . Then the molecular rates of lead and calcium from X to Y will be of the form: where A XY is a constant coefficient representing the maximum rate out of compartment X and into compartment Y, which encompasses number of cross-membrane transporters between compartments, and the average time taken by one molecule to traverse a transporter.

Absorption into blood. Gastrointestinal absorption of lead is influenced by dietary and nutritional status.
While there are a few factors which have been shown to influence lead uptake (such as food intake, iron 78 , phosphates, zinc, high fat intake, proteins, various vitamins), in this study we will concentrate on the relationship between lead and calcium competitive absorption dynamics.
An inverse relationship has been historically observed between calcium intake and blood lead concentration 79 , hence competition for a common transport protein was proposed as a potential mechanism for the lead-blood interaction at this level 69,80 . It was proposed that saturable transport mechanisms for lead may exist within the mucosal and serosal membranes and within the intestinal epithelial cell, thus affecting both intestinal absorption of lead from dietary sources, as well as blood absorption from the digestive system. These mechanisms are thought to be implemented via membrane carriers (e.g., Ca 2+ -Mg 2+ -ATPase, Ca 2+ /Na + exchange, DMT1) or facilitated diffusion pathways (e.g., Ca 2+ channel) and intracellular binding proteins for Ca 2+81, 82 . In addition, absorption of both lead and calcium from the gastrointestinal tract is enhanced by administration of cholecalciferol, which appears to involve the stimulation of the serosal transfer of lead from the epithelium, not stimulation of mucosal uptake of lead 81,[83][84][85] .
A saturable absorption mechanism supports the observed nonlinear relationships between blood lead concentration and lead intake found by a variety of studies humans and immature swine [86][87][88] , and corresponds in our model to convergence of the system to a nonlinear asymptotic steady state. To capture an overall saturable aspect to absorption, we introduced an additional gating in the model, in the form of a feedback term that slows down absorption when blood calcium is already elevated. This can also be interpreted as a protection mechanism meant to prevent an unnecessary, or even unhealthy build-up of systemic calcium. Mathematically, we considered this term to depend sigmoidally on the blood calcium levels, (with the values of = . a 0 6 and τ = 4 being kept fixed throughout the model). As discussed before, the qualitative behavior of the system is robust within a whole parameter range containing the specific values of a and τ used for our numerical experiments. Appendix A illustrates the sensitivity to varying the sigmoidal parametrs within these intervals.
Lead in blood is rapidly taken in by red blood cells, where it binds to intracellular proteins. Approximately 99% of the lead in blood is associated with red blood cells; the remaining 1% resides in blood plasma 89,90 . Studies in intact red blood cells and red blood cell ghosts suggest that there may be multiple pathways for lead transfer across the red cell membrane. While not considered the primary pathway, lead and calcium may share a permeability pathway represented by a Ca 2+ -channel 91 . Lead is extruded from the erythrocyte by an active transport pathway, likely a Ca 2+ -ATPase 92 . Altogether, our model reflects these multiple lead-calcium competitive mechanisms into nonlinear input contributions in both lead and calcium from the environmental sources p and c to the blood compartments (i.e., to the derivatives of S Pb and S Ca , respectively): ,

 
Here, the subscript ES marks that the term represents contributions from the environment (E) to the blood (S) compartment. transport in and out of soft tissue. Lead enters soft tissue from blood/serum. Mechanisms by which lead transits between blood and soft tissues have not been fully characterized 93 , but studies of mammalian small intestine suggest that lead can interact here as well with transport mechanisms for calcium and iron. Lead was shown to enter cells through voltage-gated L-type Ca 2+ channels in bovine adrenal medullary cells [94][95][96] and through store-operated Ca 2+ channels in pituitary GH3, glial C3, human embryonic kidney, and bovine brain capillary endothelial cells 97,98 . Anion exchangers may also participate in lead transport in astrocytes 93 .
To summarize these effects, we consider the rates of the two molecules from the blood compartments S Pb and S Ca into the soft tissue compartments T Pb and T Ca as The converse rates from tissue compartments T Pb and T Ca into the blood compartments S Pb and S Ca are: ( , ) ( ) can replace calcium in the calcium-phosphate salt that comprises the primary crystalline matrix of bone 99 . As a result, lead deposits are formed in bone during bone growth and remodeling and is released to the blood during the process of bone resorption 42,43 . The distribution of lead in bone reflects these mechanisms; lead tends to be more highly concentrated at bone surfaces where growth and remodeling are most active 100 , hence the bone lead distribution is age-dependent. Based on the primary calcification site, lead accumulation will occur predominantly in trabecular bone during childhood, and in both cortical and trabecular bone in adulthood 100 . Bone lead burdens in adults are slowly lost by diffusion and resorption 44,101 . The association of lead uptake and release from bone with the normal physiological processes of bone formation and resorption means that lead biokinetics is sensitive to these processes. Physiological states (e.g., pregnancy, menopause, advanced age) or disease states (e.g., osteoporosis, prolonged immobilization) that are associated with increased bone resorption will tend to promote the release of lead from bone, which, in turn, may contribute to an increase in the concentration of lead in blood [102][103][104][105][106][107] .
We model the absorption of Pb/Ca into bone as a standard nonsaturable mechanism: while for the resorption process, we will use a multiplicative parameter z (which changes with age and physiological states which have impact on bone dynamics). Resorption occurs then with release of both lead and calcium from the osseous compartment into blood, as: Effects on renal function and excretion. Granular contracted kidneys were recognized as potential effects of chronic lead exposure since the late nineteenth and early twentieth centuries. While a variety of studies documented the relationship between prolonged lead exposure and chronic nephropathy 108 , none of these studies was able to provide a conclusive proof, even though a cause and effect relationship is very likely. More recent studies have aimed to describe molecular mechanisms that may explain renal effects of lead, but with limited success. Little information is available regarding the transport of lead across the renal tubular epithelium. In Madin-Darby canine kidney cells (MDCK), lead has been shown to undergo transepithelial transport by a mechanism distinct from the anion exchanger that has been identified in red blood cells 109 . The uptake of lead into MDCK cells was both time and temperature dependent. While empirical evidence for specific transport mechanisms in the renal tubule are lacking 110 , our current knowledge suggests that (both intake and excretion) renal mechanisms are less similar to the intestinal saturable pathways of lead transfer, and more like the other soft tissue mechanisms (without the additional saturable competition) 111 . Hence the rates of the two molecules from the blood compartments S Pb and S Ca into the kidney compartments K Pb and K Ca will be also written as However, a know effect specific to the kidneys (and potentially crucial to our model and to the systemic function) is a negative feedback effect: accumulation of lead in the kidneys decreases renal function, leading to diminished excretion (of both lead and calcium, as was shown by studies relating lead poisoning to formation of kidney stones). The reabsorption rates from the kidney compartments into the blood compartments can still be considered to be simply given by but the excretion rates will incorporate the negative feedback effect as: ( , www.nature.com/scientificreports www.nature.com/scientificreports/ so that the presence of lead decreases excretion of both substances according to an exponential tail ϕ = − K e ( ) Pb kK Pb , with the parameter > k 0 describing the individual's renal sensitivity to lead toxicity. Under this scenario, lead toxicity starts gradually reducing renal function even at small levels, but can virtually shut down urine excretion when present at higher levels.

Effects on brain function and the blood-brain barrier (BBB).
There is strong evidence that adverse neurobehavioral outcomes, such as reduced IQ and academic deficits, occur at levels below 10 μg/dl. Early childhood studies of cohorts with very low exposure found a statistically significant decrease in intelligence test scores as lead levels increased from 1 to 10 μg/dl. The dose-effect relation appears to be nonlinear, with the effects of lead proportionally greater at concentrations below 10 μg/dl than above that value 7,10 , although no mechanism has been identified.
A collaboration of vascular, immune, metabolic and neural components acts as a physiological barrier (the blood brain barrier) controlling the movement of ions, molecules, and cells between the blood and the brain 112 . This protection of the neural tissue from toxins and pathogens promotes normal neural function. Dysfunction of the blood brain barrier (BBB) is associated with neurological and neuropsychiatric symptoms. Lead can penetrate the blood brain barrier and affect key processes of neural development and function, such as cell migration and synapse formation, as well as function of glial cells (which are in fact involved in BBB function). This may lead to improper brain connectivity and altered brain functions, but also further damage of the BBB, thus increasing lead transport to the brain even further, and accentuating a self-enforcing feedback loop.
Because lead and calcium compete on cross-membrane transporters, and since calcium is involved as a cofactor in many cellular processes, it is not surprising that many cell-signaling pathways are affected by lead. Lead affects virtually all neurotransmitter systems, but most mechanistic information is available on the glutamatergic, dopaminergic, and cholinergic systems, as summarized below for completion 113 .
Lead affects long term potentiation (the neurophysiological substrate for learning and storing information) in three-fold way: by increasing its threshold, by reducing its magnitude and by shortening its duration. Studies have shown that the effects of lead vary as a function of the developmental exposure period and that lead exposure early in life is critical for production of impaired LTP in adult animals. This is probably due to its action on the glutamatergic system, which has been studied both pre-and post-synaptically. Lead reduces presynaptic glutamate release, effect which is likely due to lead-related decrements in the calcium-dependent component. Some studies reported a lead-induced postsynaptic increase in number/density of glutamate receptors, but results regarding the effects of lead on postsynaptic glutamatergic function have generally been inconsistent.
Studies in animals also report effects of lead on nigrostriatal and mesolimbic dopamine systems regarding receptor binding, dopamine synthesis, turnover, and uptake (although the effects of these on cognitive function have been inconclusive). Exposure to lead induces numerous changes in cholinergic systemic function, which also plays a role in learning and memory processes. It was shown that lead blocks evoked release of acetylcholine and diminishes cholinergic function, in both central and peripheral synapses. This mechanism may involve calcium dynamics as well, in that lead reduces acetylcholine release by blocking calcium entry into the terminal and prevents sequestration of intracellular calcium by organelles, which results in increased spontaneous release of the neurotransmitter.
One pathway that has been studied in more detail is the activation of protein kinase C (PKC). PKC is a serine/ threonine protein kinase involved in many processes important for synaptic transmission such as the synthesis of neurotransmitters, ligand-receptor interactions, conductance of ionic channels, and dendritic branching. One of several calcium-dependent forms of PKC is a likely target for lead neurotoxicity; studies in vitro showed that it is neuron-specific and is involved in long-term potentiation, spatial learning, and memory processes. However, studies in rats exposed to low lead levels have shown few significant changes in PKC activity or expression, suggesting that the whole animal may be able to compensate for lead PKC-mediated effects compared to a system in vitro. PKC induces regulation of an astrocytic gene (GFAP). Astrocytes along with endothelial cells make up the BBB. Studies in rats exposed chronically to low lead levels have reported alterations in the normal pattern of GFAP gene expression in the brain, and the most marked long-lasting effects occurred when the rats were exposed during the developmental period. It appears that premature activation of PKC by lead may impair brain microvascular formation and function, and at high levels of lead exposure, may account for gross defects in the blood-brain barrier that contribute to acute lead encephalopathy. The BBB normally excludes plasma proteins and many organic molecules, and limits the passage of ions. With disruption of this barrier, molecules, ions and water enter the brain more freely and, given the slow lymphatic drainage, lead to edema and increased intracranial pressure. The particular vulnerability of the fetus and infant to the neurotoxicity of lead may be due in part to immaturity of the blood-brain barrier and to the lack of the high-affinity lead-binding protein in astroglia, which sequester lead 48 .
Although the precise mechanism for the inhibition of glutamate release by lead is not known, it is consistent with lead preventing maximal activation of PKC, rather than lead blocking calcium influx into the presynaptic terminal through voltage-gated calcium channels. Moreover, glutamate release was shown to have a U-shaped response: it was inhibited in rats treated with the lower lead doses, but not in those exposed to the higher concentrations of lead. Although speculative, this was interpreted as lead at the higher doses mimicking calcium in promoting transmitter release and overriding the inhibitory effects of lead that occur at lower lead levels 48 .
Our model aims to reproduce these complex effects in a simplified way: we assume that, even though lead tries to mimic calcium by using the same molecular mechanisms, the existence of multiple molecular paths for calcium through the healthy BBB confers it the ability to "differentiate" to some extent between lead and calcium molecules. This ability decreases with the accumulation of lead in the brain. We incorporate these effects of the BBB into our rates from the blood to the brain compartments as: www.nature.com/scientificreports www.nature.com/scientificreports/ Here, ξ is a decreasing function of the brain level of lead B Pb , which we considered for specificity to be ξ = − − Be 1 sB Pb . Here, the parameter B represents the BBB filtering efficiency, and s represents the sensitivity of the BBB to the neurotoxic effects of lead. Briefly speaking, an elevated brain level of lead B Pb increases ξ towards one, and subsequently not only increases the probability for lead to traverse through the Pb/Ca transporters, but also shifts the effective interval of permeability to both. The significance of both parameters to the dynamics is further interpreted in the Results section.
The reversal rates from the brain to blood compartments are not gated by the BBB, and will be considered to simply follow the non-saturable mechanism described by:

Lead-calcium model
Combining all these coupled compartments, our model will have the following form (where the variables of the functions, as described in the paragraphs above, were omitted here for clarity: where the functional dependences are defined in the previous section, based on existing qualitative information on the corresponding molecular mechanisms. Before starting the model analysis, let us establish a unit convention for our measurable quantities. In our simulations, time is expressed in weeks. In line with the typical units used in empirical studies, the variables, representing concentrations of lead or calcium, are all expressed in mg/kg, for consistency of the equations. For the fluid compartment S, 1 mg/kg is virtually equivalent to 1 mg/l (if we approximate the volumetric density of blood with that of water). The same applies to the input parameters p and c, if we consider the source of lead/calcium intake to be a liquid. The units for the sigmoidal parameters are such that the nonlineaity  is nondimensional, and the flow rates A XY are measured in mg/kg/weeks. In this study, we will consider, for simplicity, that all relevant parameters A XY are identical. These may have in reality different values for different compartment pairs (as specified by the lower indices). Since we do not yet have quantitative empirical estimates of the permeability parameters, and since in this project we do not explore the mechanisms and effects of increasing these maximal rates, we will normalize all these coefficients to one for the rest of the paper, in order to fix our ideas and simplify the analysis.

Results
We first investigate the dependence of the system's long term behavior on the lead versus calcium intake alone. We will then study how this dependence is affected by perturbations in other system parameters, in particular age (as reflected in the bone resorption parameter z, as well as BBB efficiency B and sensitivity s), and kidney sensitivity k to lead toxicity.
Dependence on lead intake p. Throughout this section, the calcium intake is kept fixed ( = c 10), as well as all other system parameters (as specified in the figure captions). The parameter p (lead intake) was allowed to increase from small (close to zero) to large. The overall calcium intake was specifically fixed to a large enough value so that c remains larger than p throughout this analysis, as common sense and biology would suggest.
The overall effect observed in all compartments ( = X S T B K O , , , , ) was that the lead content X Pb generally increases, and the calcium X Ca content generally decreases with increasing lead intake p. This is not surprising, given the competitive lead/calcium mechanism on which our model relies. This is consistent with empirical With the unit conversion of mg/m 3 = mg/l, the behavior also matches data and simulations by O'Flaherty in subjects exposed to long-term air contamination, as shown in Williams et al. 116 for a 30 year old male exposed to air lead for 10 years, as well as for a group of 30 subjects with long term exposure 34,117 (also see 43 , Figs 7 and 8 in the reference).
Interestingly, however, the dependence on p is not just a simple monotonic trend: the equilibrium curve with respect to p exhibits two saddle node bifurcations. To fix our ideas, we plot four different projections of this equilibrium curve in Fig. 3: the B Pb and B Ca components (in gray and cyan, respectively) in Fig. 3a, and the S Pb and S Ca components (in gray and cyan, respectively) in Fig. 3b. Figure 4 further illustrates S S ( , ) Pb Ca and B B ( , )

Pb
Ca phase space slices for representative values of p, as we will further discuss in the next paragraph.
For very low lead intake p, the system has a unique locally stable equilibrium. Slightly increasing p gradually increases the lead equilibrium content in all compartments, and lowers the calcium equilibrium content in all compartments. However, at = .
p 0 28, a second locally stable equilibrium appears, via a saddle node bifurcation (shown as a green square along the equilibrium curve). We will comment more on this bistability window in the next paragraph and associated figures. Bistability ends at = .
p 0 85, via a second saddle node bifurcation (shown as a yellow square along the equilibrium curve). For values of p larger than this bifurcation value, the system again has a unique stable equilibrium, quantitatively different, however, than before bistability occurred.
Interestingly, the two bistable equilibria differ only in their brain components B Pb and B Ca , as shown in Fig. 3a. Figure 3b illustrates the situation in the blood components, where S Pb and S Ca are identical between these equilibria (this also being the situation more generally, for all other compartments except for the brain). Bistability allows the system to converge to either stable equilibrium, based on its initial conditions. To fix this idea, the panels of Fig. 4b illustrates S S ( , ) Pb Ca and B B ( , ) Pb Ca phase space slices for values of p before, within and after the bistability window. To investigate how the long-term behavior of the system depends on prior short-term exposure to lead and ingestion of calcium, we varied the blood components of the initial conditions (which act as short-term storage), and fixed the other initial components to the low value 0.1. For the initial values of S S ( , ) Pb Ca , we considered a grid of size 1 for the 2-dimensional square × [0, 10] [0, 10]. The phase slices show that, for the bistability parameter, some of the solutions converge to the high B Pb and low B Ca equilibrium, and others to the low B Pb and high B Ca equilibrium. It is worth noticing again that (1) bistability can only be detected when looking at the brain components of the system; (2) when in the bistability window, whether the system converged to one versus the other equilibrium can only be distinguished by looking at the brain components.
To better illustrate how localized changes in the initial conditions may prompt the system to converge to either equilibrium when within the bistability window, we show in Fig. 5 all components of the solutions for a bistable value of p, for two sets of initial conditions. All solutions, after a transient phase reflecting the system's kinetics, converge asymptotically to an equilibrium, as expected from the bifurcation diagram. In all panels, the time units are weeks, and the steady state blood concentrations of both lead and calcium are measured, as mentioned in the Methods section, in mg/l. The top panels consist of temporal evolutions for all variables, when initiated at a low value (0.1), and the bottom panels show the same temporal evolutions when the initial blood lead only was elevated to = S 4 Pb . An initial spike in blood lead (when p is in the bistability bracket) affects, counter-intuitively, not the long term blood or bone levels of lead, but the long-term brain levels of lead and calcium (raising the brain lead, and diminishing the calcium brain levels). Notice that, in both scenarios, the system is settled towards a steady state of around 0.3 mg/l = 30 μg/dl within the first two to three weeks, which is in agreement with empirical evidence from a 25 year old man receiving daily lead intake over a prolonged period, as documented in Cools et www.nature.com/scientificreports www.nature.com/scientificreports/ al. 118 , and also modeled by O'Flaherty 43 . Notice that one could not distinguish between the two different steady states of the system by looking only at the kinetics of the blood compartment.
Moreover, empirical data from both single subjects 119 and groups of subjects 120 show a monotonic, close to diagonal dependence of kidney lead long term concentration to blood lead concentration, which is in line with our plot (high kidney lead occurs in conjunction with equally high blood lead and low kidney lead occurs in conjunction with equally low blood lead). Data suggest that this is not the case in the blood-bone lead concentration dependence 121 , which is confirmed by our graphs. More qualitative effects will be further discussed in the next section.
Dependence on calcium intake c. Figure 6 shows that, as one would expect, the equilibrium content at of blood calcium increases with calcium consumption c, and that the equilibrium content of blood lead decreases with calcium consumption. This trend does not depend on the lead intake p, appearing qualitatively consistent between the three different lead exposures we considered: low ( = .
p 0 1, blue curve), medium ( = . p 0 5, magenta curve) and high ( = . p 1 5, brown curve). Although the curves for elevated lead intake ( = . p 0 5 and for = . p 1 5) exhibit saddle node bifurcations (leading to sizable bistability windows for c), the blood components of the two coexisting equilibria are identical (i.e., the two equilibria could not be distinguished from each other if one only compared their S Pb and S Ca components of the steady state. This is in fact the situation for all other system compartments (not shown in the figure), expect for the brain, which represents the only compartment that can distinguish between the two coexisting locally stable equilibria.
For the brain components, the dependence of brain calcium/lead dynamics on calcium intake is more complicated. This dependence is simply monotonic for a fixed small lead intake ( = .
p 0 1, blue curve in both panels a and b): brain lead levels are small, and decrease with c; brain calcium levels are high, and increase with c. For high levels of lead intake ( = .
p 1 5, brown curves), brain calcium levels are small, but slightly increasing with c for a wide range of c. At a critical value for c (saddle node point ∼ c 18, marked with a green square along the brown curve), the system gains access (depending on the state of the system) to a second locally stable equilibrium branch, with much higher calcium, increasing even further with c. In turn, brain lead levels are high for high www.nature.com/scientificreports www.nature.com/scientificreports/ levels of lead ( = . p 1 5). While they raise significantly with c at low c values, increasing c more severely will eventually start diminishing them. The saddle node value ∼ c 18, brings access to the second stable equilibrium branch, with small B Pb decreasing even further with c.
For intermediate values of p (pink curve), the brain calcium is close to zero for no calcium intake c, as expected. Then the curve exhibits a bump at relatively low values of c, followed by a slight dip, and a recovery (lower stable branch of the pink curve in Fig. 6b). Hence brain calcium starts performing well with increasing c, even in the presence of moderately high p, but unfortunately this also reinforces the brain lead component. The maximum point of this increasing trend occurs close to the first saddle node bifurcation value at ∼ c 4. Past this point, within the bistability window for c, the system may converge, based on its history (i.e., initial conditions) to one of two stable regimes: either a productive one, of high B Ca and low B Pb (with B Ca further increasing and B Pb further decreasing with c) or a detrimental one, of low B Ca and high B Pb (with B Ca further decreasing and B Pb further increasing with c). When c is increased even more, and the system leaves the bistability window, only the productive equilibrium remains, although the lead/calcium balance will not improve significantly when further increasing c.
The bistability window appears to be a crucial and unexpected feature that breaks the intuitive monotonicity result. To study the onset of this window in terms of both lead and calcium intake, we followed the saddle node curve in the parameter plane p c ( , ), shown in Fig. 7. The figure illustrates a subset of the 2-dimensional parameter slice p c ( , ) (with the other parameters fixed to the same values as those specified in the previous figure captions). The two branches of the saddle node curve, meeting at a cusp point, delimitate the region in which the system has bistability (area "inside" the curve) versus a unique stable equilibrium (area "outside" of the curve). Read from left to right, this plot shows that, for fixed calcium intake, increasing lead intake will always bring the system through a bistability window (except for very low calcium c, where the transition from low to high calcium/lead balance occurs without phase transitions). Increasing dietary calcium c for a given level of lead contamination may not necessarily amount to higher calcium and lower lead content in the system compartments (as shown in Fig. 6), but rather will increase the size of the window in which the long-term calcium/lead distribution in the system depends crucially on its initial conditions (history of exposure). This is important, and will be further discussed in Section 3.
In the next subsections, we investigate the effect that other physiological parameters (z, k, B and s) have on altering the details of the system's balance and long-term dynamics. We then further interpret and contextualize these details in the Discussion section.
Dependence on kidney sensitivity k. One of the significant effects of lead toxicity outlined in the Introduction and Methods sections is that on kidney function. In our model, we represented this effect through a parameter k, so that the efficiency of kidney-filtered excretion is modulated by the amount of lead accumulated in   Fig. 4): low lead intake ( = . p 0 1, blue); medium lead intake, corresponding to p within the bistability window ( = .
p 0 5, magenta); high lead intake ( = . p 1 5, brown). Saddle node bifurcations are marked as squares along these equilibrium curves. Locally stable branches as shown as solid lines, and unslable branches as dotted lines. The other system parameters are fixed, as before, to: = z 1, = .  produced by potentially different tuning of physiological mechanisms of kidney vulnerability, but which may also be due to unrecognized sources of exposure. The author proposes that a physiological model may help with understanding which is the case, and generate more reliable predictions. Below, we analyze potential effects on dynamics produced by varying k in our own model. Figure 8 shows a quantitative difference between the case of low and that of high kidney sensitivity, in that the steady state of the system is affected as the lead intake p is increased. Increasing k has a very intuitive, straight-forward effect on the kidney compartments K Pb and K Ca , raising dramatically the levels of both renal lead and calcium when lead intake p is increased (bottom panels). Notice that the severe renal functional impediment introduced by raising k tenfold transcends the normal trend of K Ca decaying with higher p, and instead leads to a large renal storage of calcium, as well as lead. Increasing k has a less intuitive effect on the brain lead/calcium dynamics. While higher sensitivity k leads to substantially increasing the accumulation of brain lead B Pb when increasing the intake p, the level of brain calcium remains unaffected, even with wide variations of k. This will be further interpreted and connected to potential clinical aspects in the Discussion.
The effects of varying calcium intake c are relatively robust to increasing the kidney sensitivity k. Even with a ten fold difference in k, the brain components of the steady states are almost identical. As one would expect, increasing k slightly increases calcium retention in the kidneys (Fig. 9, bottom right panel). However, the truly sizable difference can be observed in kidney lead retention in the low calcium regime, when continuing to decrease calcium past the saddle node threshold allows lead to substantially build up in the K Pb compartment (Fig. 9, bottom left panel).
Notice that, in general, while there are substantial quantitative differences in specific components of the system's steady state in response to varying k, the bifurcation points are surprisingly robust to these changes. This effect can be observed for the examples in Figs 8 and 9 (the saddle node bifurcations occur at approximately the same values for p and c when k is changed), but is also represented as a more complete picture in Fig. 10. k 0 05 (blue curves) and higher sensitivity = . k 0 5 (red curves). Each panel shows a different projection of the same two equilibrium curves. We chose to illustrate the kidneys components K Pb and K Ca , since they are most likely to be affected, and the brain components B Pb and B Ca , since they are the primary object of our study. The other system parameters are fixed, as before, to: = c 10, = z 1, = s 1, = . www.nature.com/scientificreports www.nature.com/scientificreports/ Notice how, as k increases, the bistability window with respect to p conserves it size, but slightly drifts to lower p (Fig. 10a), while the bistability window with respect to c slightly drifts to higher c (Fig. 10a).
The same effect can be visualized in the p c ( , ) parameter plane. Figure 11 shows in blue the two branches of the saddle node bifurcation curve obtained for = .
k 0 05 (as in Fig. 7), colliding at a cusp point (marked in green). It also shows in red the same two branches for = .
k 0 5, colliding at a cusp point marked in yellow. The two curves are almost indistinguishable, illustrating the robustness of the bistability region with respect to k.
Dependence on bone resorption z. The bone resorption rate (as a marker of age or of other physiological states such as pregnancy or osteoporosis) is dependent in out model on the parameter z, which alters the sensitivity to accumulation of calcium in the osseous tissue: increasing z increases the response and prompts faster resorption at lower bone calcium content. Below, we will analyze the evolution of the system's steady state for different degrees of calcium and lead intake, under different regimes for z. Figure 12 shows the dependence of the osseous equilibrium components on lead intake p, for two different levels of bone resorption z. As expected in light of the previous section, the bone level of lead generally increases with lead intake, and the bone level of calcium generally decreases with lead intake. One can also notice, however, that the the equilibrium levels of both bone lead and calcium are more pronounced for lower levels of bone resorption, allowing both to accumulate faster in the osseous compartment ( = z 1, red curve), and less pronounced with higher resorption, as calcium (and lead together with it) leaves the bones at a faster rate, transitions through other compartments and is eventually eliminated.
Bistability occurs for an interval of p which appears independent on z. The bone components of the coexisting locally stable equilibria are identical, making bistability irrelevant in the bone compartment subspace. As before, bistability is quantitatively significant for the brain projection of the equilibrium, which evolves as shown in Fig. 3 (i.e., exhibits a high lead, low calcium state and a low lead, high calcium state). The brain levels of lead and calcium, however, do not change with z, and the equilibrium curve is identical for all resorption levels. www.nature.com/scientificreports www.nature.com/scientificreports/ The behavior is equally straight forward when considering changes in response to increasing the calcium intake c. As shown in Fig. 13, the bone lead content decreases asymptotically to a value close to zero, and the bone calcium content increases from zero to an asymptotic value, when the calcium intake c is increased. As expected, the bone accumulation is more pronounced in both lead and calcium for smaller values of z (lower resorption rate). One can notice the large bistability window with respect to c, with two coexisting locally stable equilibria with all identical components except for the brain compartments, illustrated in the bottom panels of Fig. 13 (one characterized by high lead and low calcium, and the other by low lead and high calcium). The onset an offset of the bistability window do not depend on the parameter z. In the Discussion section, we further interpret the significance of this result in the context of age effects on the long-term behavior of the system.
Dependence on Blood Brain Barrier (BBB) efficiency B and sensitivity s. In this section, we investigate the effects of the BBB filter on the brain Pb/Ca balance. To do this, we consider two key parameters. The first is the efficiency ≤ ≤ B 0 1 of the BBB in discriminating between "good" materials, that promote brain function (e.g., molecules such as calcium, which should be allowed through the filter) and "bad" components, that are neurotoxic (e.g., lead molecules, which should be gated). The second parameter s represents the BBB sensitivity to lead neurotoxicity (that is, the dependence of the BBB accuracy on the brain's existing neurotoxicity). When , but as B Pb builds up, the exponential term decreases, gradually weakening the Pb/Ca gate away from its natural efficiency described by B. In absence of brain lead, an ideal filter ( = B 1) would completely filter out lead, and only allow calcium to pass through the BBB. In our simulations corresponding to  Figure 11. Saddle node bifurcation curves in the p c ( , ) parameter plane, for two different kidney toxicity levels = . k 0 05 (blue curve) and = .
k 0 5 (red curve). The cusp points are marked along the curves as a green and respectivelly yellow diamond. For each k, the bifurcation branches enclose the p c ( , ) region of bistability, which appear thus robust with respect to k. The other system parameters are fixed, as before, to: = z 1, = s 1, = . www.nature.com/scientificreports www.nature.com/scientificreports/ healthy adult systems, we work with values of B close to the optimal value, illustrating a fully developed and efficiently tuned BBB system. However, allowing the sensitivity s to increase effectively acts as lowering B, and the efficiency drops away from the optimal tuning.
We focus on illustrating the behavior of the brain components B Pb and B Ca , which are the components primarily affected by changing B and s. Figures 14 and 15 show the effect on brain dynamics of increasing lead and respectively calcium intake, for two different BBB efficiency regimes: high efficiency ( = . B 0 9, equilibrium curves shown in red) and low efficiency ( = . B 0 7, curves shown in blue). Broadly speaking, decreasing the efficiency B leads to increasing the brain lead steady state for a wide range of lead intake p (Fig. 14a) and calcium intake c (Fig. 15a). Figures 14b and 15b show that decreasing the efficiency B leads to lower brain calcium steady states, throughout both the lead and the calcium intake range.
However, the effect of varying B is a little more complicated than that. A more accurate description of this effect is to phrase it in terms of a dynamics shift. Decreasing B produces approximately a translation in the equilibrium curves, which results in a shift of the saddle node bifurcations and of the bistability window: to the left, in Fig. 14, and to the right, in Fig. 15.
We represent this shift in behavior with respect to p and c simultaneously in Fig. 16, by tracking the change in the position of the saddle node curve delimiting the bistability window in the p c ( , ) parameter plane, as the filter efficiency changes from = .
B 0 9 (red curve) to = . B 0 7 (blue curve). Notice that for higher efficiency B, one needs to additionally increase p and/or decrease c in order to create the potential for bistability.
This relationship between B and the position and size of the bistability locus has subtle, but important effects on the model dynamics. Our plots essentially imply that, for a low efficiency filter B, lower doses of lead intake p have the effects that higher doses would have when applying a more efficient filter. The effect is even more sizable when varying calcium intake c: higher doses of c produce effectively to the same steady state levels for which lower doses of calcium would be sufficient if operating with a more efficient filter. This implies, for example, that for higher efficiency the low B Pb and high B Ca equilibirum becomes reachable at lower values of c, and that the high B Pb and low B Ca steady state disappears at lower c values, both of which are potentially beneficial effects when considering addressing lead toxicity by increasing calcium intake. www.nature.com/scientificreports www.nature.com/scientificreports/ Varying the sensitivity s of the BBB to lead neurotoxicity effectively reflects into varying the efficiency B. We discuss this dependence separately for three reasons. First, the dependence of BBB "efficiency" on s is nontrivial, and we expect the results to reflect this situation. Second, the sensitivity s acts in conjunction with the lead neurotoxic levels, therefore introduces a feedback which represents one of the weak points that lead pharmacokinetics  www.nature.com/scientificreports www.nature.com/scientificreports/ triggers in its path through the system. Third, we consider efficiency and sensitivity to reflect two different mechanisms embedded in the BBB, which may act together with combined or opposite effects. A low efficiency system, for example, may perform better than a more efficient one under lead toxicity conditions, due to an increased resilience (low sensitivity) to neurotocixity. Below, we discuss the perturbations in the system behavior produced by varying s. Figure 17 illustrates the steady state levels of both brain lead and calcium. Unlike the effect of decreasing the efficiency B, which only produced a slight shift in these equilibrium curves, the effect of varying s is more dramatic. As expected, higher sensitivity generates higher lead levels B Pb and lower calcium levels B Ca . However, this trend is quantitatively uneven along the lead intake range, as well as within the range considered for s. First, notice that prominent effects of varying s occur in the lower sensitivity regime, and they saturate as s increases: the switch between the blue and the red curves, when changing s from low ( = . s 0 5) to medium ( = s 1), is visibly larger than that between the red and the green curves, when changing s from medium to large = .
s 1 5). Secondly, notice that the effect of varying s is minor at large lead intakes, but substantial in the medium and low p range, where bistability occurs.
If the BBB has unusually low sensitivity then the steady state levels of B Pb increase monotonically with p, and the steady state levels of B Ca decrease monotonically with p (blue curves). This effect is not surprising in a system where small increases in neurotoxicity leave BBB function relatively undisturbed. The dependence on p is also monotonic when the BBB sensitivity is unusually high, although it is quantitatively very different: the steady state B Pb is already high for low p and only slightly increases with p; and the B Ca steady state is extremely small and virtually constant with respect to p. This is not unexpected either, since a in a system with a filter highly sensitive to neurotoxicity, even small p levels lead to effects on brain dynamics comparable to those inflicted by more substantial lead intakes (that is, high long-term B Pb and virtually no B Ca ).  www.nature.com/scientificreports www.nature.com/scientificreports/ It is for the intermediate values of s that the dependence of lead intake is the most complex, where simple monotonicity of equilibria is replaced by the emergence of saddle node bifurcations, and the subsequent bistability window. Figure 18 illustrates the dependence of equilibria on calicium intake c, as the sensitivity s is varied. At low sensitivity = .
s 0 5 (blue curves), both B Pb and B Ca components of the steady state are monotonic, following the intuitive trend: a higher calcium intake lowers the equilibrium brain lead, and increases the equilibrium brain calcium. When increasing s, the bistability window forms and evolves with increasingly larger values of s, as follows.
For medium and high BBB sensitivity values s, the equilibrium curve -illustrated in red ( = s 1) and respectively in green ( = . s 1 5) in Fig. 18 -exhibits saddle node bifurcations and bistability with respect to varying c. Within the bistability window, the system has access to two distinct steady states, one characterized by high brain lead and low brain calcium, the other corresponding to low brain lead and high brain calcium. Small changes in initial conditions (e.g., a slightly higher exposure to lead toxicity in the past) may swap between attraction basins, and prompt the system to evolve toward the ''bad" instead of the "good" long-term prognosis.
When pushing c below the lower saddle-node entry to bistability, the low lead/high calcium equilibrium is lost, and only the high lead/low calcium equilibrium survives. This branch exhibits the behavior described in Section 2.2, suggesting that calcium control may exercise counter-intuitive effects in the low calcium range: both B Pb and B Ca steady state values decrease with lowering c, except for a slight bump in brain calcium in the low c range, before B Ca eventually decays to zero for zero calcium intake. This "inverted" effect (both the size of the calcium bump and the slope of the lead decay) attenuate with higher sensitivity (as shown by the differences between the red and the green curves, for medium and respectively high BBB sensitivity). When pushing c past  www.nature.com/scientificreports www.nature.com/scientificreports/ the upper saddle-node exit from bistability, the high lead/low calcium equilibrium is lost, and only the low lead/ high calcium one persists, suggesting that large enough calcium intakes can counteract the effects of moderate lead toxicity (all graphs correspond to = .
p 0 5 lead intake). Increasing sensitivity extends the bistability window, moving the exit saddle node to significantly higher values of c (e.g., for = . s 1 5, this value was too high to represent within the normal domain for c). It is easy to understand why higher BBB sensitivity can lead to such an effect: the potential for a "bad" steady state remains available to the system, for a basin of vulnerable initial conditions, even under high calcium treatment (treatment which would efficiently eliminate this potential for a system with lower sensitivity). For each sensitivity level s, the critical window of lead intake p for which this bistable behavior occurs is illustrated in Fig. 19a. For each sensitivity level s, the onset and offset of bistability with respect to c are illustrated in Fig. 19b. Figure 20 shows the simultaneous effects that BBB sensitivity has on the p c ( , ) bistability locus.

Discussion
Specific comments on the model. In our model, we introduced basic molecular assumptions of competitive transport of lead and calcium, to describe the compartmental dynamics of these two elements in humans.
Our primary interest was to observe changes in dynamic regimes in response to varying the doses of both lead and calcium intake, and to further analyze how these transitions are altered when introducing constraints on different system parameter related to age, physiological and disease states (such as bone resorption rate, kidney sensitivity to lead, blood brain barrier efficiency and sensitivity to neurotoxicity). This is a phenomenological model, aiming to prove that it is possible in principle to use a nonlinear compartmental model to simulate and  www.nature.com/scientificreports www.nature.com/scientificreports/ better understand the effects of lead-induced toxicity. Some of our results are straight forward, supporting existing empirical evidence of the dose-effect relationships. Other results are more complex, and suggest that further investigating the competitive pharmacokinetics between lead and calcium (and perhaps more generally, between lead and other interacting molecules) may represent a promising avenue of researching lead neurotoxicity.
We first studied the effects of varying the input the the system. Higher lead intake resulted in simply increasing the lead content and decreasing the calcium content, in all compartments except for the brain. In the brain, increasing lead intake may initially result in increasing the potential for a higher lead and lower calcium steady state, by entering a bistability window (with two potential asymptotic outcomes, a "good" one, with low brain lead and high brain calcium, and a "bad" one, with high lead and low calcium). If the lead intake is increased even further, the system exits the bistability, and only the high lead, low calcium steady state persists. Crossing the bistability window represents a relatively mild transition between regimes, which does not suddenly degrade the system performance, and allows it to keep performing well under certain circumstances (e.g., history or prior exposure), even under rather significant increments in p. The phenomenon emphasizes the importance of the system's history of toxicity when attempting to predict its response to newer doses. A similar sequence of phase transitions appears in response to increasing calcium intake. We showed that, irrespectively of the lead ingestion p, an intense enough treatment of the system with calcium will always eventually reduce brain lead and elevate brain calcium. However, when p is high, the necessary calcium dose may not be physiologically plausible. Interestingly, our model also suggests that when exposure to lead is very low, administering too little calcium may in fact be more detrimental to the system (enhancing brain lead) than no calcium intake.
One potentially important observation around these results is that the effects of lead and calcium ingestion on the brain can only be observed in the brain, and cannot be measured or inferred from their levels or from their ratio within any of the other compartments. A second important observation is that, while the calcium/lead transport mechanism is competitive between all compartments, and was incorporated as such in our model, this does not mean that low lead in a compartment is necessarily associated with high calcium in that compartment, or the other way around. The model includes gating and feedback mechanisms, which lead to counter-intuitive dependences of the long-term behavior of the system with respect to the doses of lead and calcium ingested, especially when the system's history (initial conditions) comes into play.
We then investigated how variations in excretion dynamics may affect the systemic asymptotic levels of lead and calcium. We focused on modeling and understanding the negative feedback that lead toxicity induces by acting as a blocker of renal function; this prevents normal lead excretion and contributes to higher accumulation of lead in the system, and hence to higher toxicity. As expected, the renal functional impediment introduced by raising k facilitates a large renal storage of both calcium and lead, especially in response to increasing the intake of either, and a sizable built-up of renal lead toxicity even in response to low lead doses (augmented when the system functions with very low calcium intake c). Higher sensitivity k also lead to a pronounced accumulation of brain lead B Pb when increasing the intake p; however, quite surprisingly, the levels of B Ca in response to increasing p remained generally unaffected, even with wide variations of k. Overall, while the kidney sensitivity gating may explain well known effects related to renal function (such as calcium storage potentially leading to kidney stones, or renal lead storage that may contribute to further enhancing kidney damage), and may even contribute to tuning brain toxicity, it does not seem to be responsible for controlling brain calcium levels, which remained robust in our model for a large range of sensitivity values.
Since one primary goal of our modeling project was to understand the effect of lead on early brain function, we next investigated the potential effects of age on the system dynamics, in particular on the lead/calcium balance. The age-related aspects built into the model design are the bone resorption parameter z, and the blood brain barrier parameters B and s (the first representing the intrinsic BBB efficiency, and the second -the dependence of its efficiency on brain toxicity).
Our analysis of equilibiria and bifurcations with respect to z suggests that, as one would expect, variations in bone resorption rates affect primarily the osseous compartments. Higher bone resorption (typically associated with older age or physiological states like pregnancy, osteoporosis or menopause) readily lead to loss of both bone calcium and bone lead. In a healthy system, the resulting excess of these in the blood compartment rapidly transitions to the kidneys and is excreted. In a system with poor renal regulation, the diffusion of bone storage into other compartments may contribute to increasing active toxicity in the system, an avenue which we have not particularly investigated in our model analysis. When assessing the effects on the brain compartment, we noticed that variations in z did not produce in our model any direct impact on lead/calcium dynamics, leaving the long-term levels of brain lead and calcium unaltered, as well as the points where phase transitions occurred, implying that the bistability windows were robust with changes in bone resorption regimes. Hence we can overall conclude that the model predicts that neural deficits do not occur as a direct consequence of increased bone resorption brought on by age or physiological states.
The BBB was represented in our model by two parameters: B, describing the efficiency of the BBB filter in a regular, nontoxic environment, and s, describing the sensitivity, or vulnerability of the BBB function to brain lead neurotoxicity. The BBB is a brain component that continues to develop throughout childhood. While there is a lot that is not known with respect to which aspects of the development occur in what order, a safe modeling assumption is that the efficiency of the BBB filter increases rapidly from lower to higher values during infancy and childhood, and that the sensitivity decreases over this age span.
We first noted that, unsurprisingly, stronger BBB efficiency is an advantage for the system, limiting the buildup of lead and allowing higher storage of calcium in the brain compartment. However, a more important, qualitative observation relates to the change in position of the bistability window as B is increased. We observed that, for lower efficiency B (which can be interpreted as an immature BBB): (1) lower doses of lead intake p have the effects that higher doses would have when applying a more efficient filter and (2)  www.nature.com/scientificreports www.nature.com/scientificreports/ efficient filter. In particular, the on and offset of bistability are affected. For a system functioning in a regime with a sustainable low B Pb /high B Ca equilibrium, a higher lead dose is required for a more efficient BBB than for a less efficient BBB to transition to the bistability regime, where the system gains access to the unsustainable high B Pb / low B Ca equilibrium. Down the line, a higher additional lead dose is required to transition our of the bistability window, so that the system loses access to the sustainable equilibrium, independently on the initial conditions. These effects on dynamics are more pronounced when varying the BBB sensitivity s to lead neurotoxicity. The sensitivity s can be assumed to be generally decreasing with age, with higher s values representing an immature or deficient BBB filter (i.e., young age or physiological conditions affecting BBB integrity). We showed that the sensitivity s also affects the size and position of the bistability window. In our simulations of a typical resilient system, with low sensitivity to lead neurotoxic effects, bistability does not occur, and the dependence on the lead and calcium intake is simply monotone. Bistability emerges at a critical sensitivity level, and increasing sensitivity further will significantly affect the onset and offset of bistability in terms of the lead/calcium intakes. For a chronological, "age progression, " one can notice that for an immature (young, very sensitive) BBB,the first saddle node bifurcation occurs at lower lead intakes (opening accessibility to the "neurotoxic, " high lead/low calcium steady state) and so does the second saddle node bifurcation (closing accessibility to the productive, low lead/high calcium steady state). As the BBB matures, sensitivity decreases and bistability narrows and moves to increasingly higher lead intakes, until it may even disappear, if the BBB desensitizes (with the inherent advantages and disadvantages). For an immature BBB, bistability onset also and offset also occur at higher calcium intakes, suggesting that a more substantial calcium treatment is necessary to address neurotoxic effects than in younger subjects, and that the potential for the neurotoxic steady state persists in a sensitive system even when aggressively treated with high doses of calcium. This situation is resolved when the BBB matures, and the system moves in a different s regime, in which the the access to the neurotoxic state can be closed up by lower calcium doses.
The existence of the bistability regime and its dependence on brain development may therefore represent the underlying basis for a potential mechanism explaining the nonlinear trends documented empirically in children with low level exposure. Low level lead exposure typically refers to blood lead at concentrations less than 30 μg/ dl, and may not produce any measurable physiological symptoms. It has been argued, however, that low level contamination can lead to significant and long term cognitive and behavioral deficits in developing children 10 . A few known studies converge to show a monotonic relationship between blood lead concentration and the size of cognitive deficits. Using data from a cross sectional study of more than 5,000 U.S. children ages 6-17 years, Lanphear et al. 77 focused on subjects within blood lead concentration under 10 μg/dl and related higher blood lead to poorer cognitive performance (measured through nonverbal reasoning, short-term memory, arithmetic, and reading skills). A cohort study of children living in Rochester, as well as an analysis of a subsample with very low exposure, also found that children's intelligence test scores declined as lifetime average blood lead levels increased from 0 to 10 μg/dl 76 . An analysis of data from 48 children from Boston 75 , with blood lead concentrations less than 10 μg/dl, mirrored these results 74 . For confirmation, data were used from 1,333 children who participated in seven international population-based cohort studies. The IQ of the participants, followed from birth or early infancy until 5-10 years of age, was again shown to decrease with increasing low level concentrations. Interestingly, all three studies also supported a nonlinear dose-effect relationship, in the sense that the effects of lead appeared to be proportionally greater at lower concentrations.
In our model, one may regard the levels of brain calcium as a rough measure of cognitive performance, since calcium is a main promoter of neurotransmitter function, which is the basis of cognition. More precisely, while there may be neural effects of lead contamination which do not necessarily involve calcium dynamics (as mentioned in Section 1.5.5), we may safely assume, for simplicity, that low brain calcium concentrations lead to poor cognitive performance, and higher concentrations are associated with improved performance. In this context, the model generally supports both the monotonicity and the nonlinearity of the dose to effect relationship observed in data (as described above). Notice first that the monotonicity is satisfied at the more basic level of the relationship between intake lead and blood lead concentration after long term exposure. However, the model predicts this dependence to be linear. The nonlinearity emerges in the relationship between dose and brain/calcium lead concentrations. These are difficult to measure directly in human subjects, hence we will use the predicted brain calcium levels as the signature of cognitive performance, in order to further compare our results with data.
Our model suggests that, for parameters corresponding to a healthy adult (with high BBB efficiency and low BBB sensitivity), there is a range of low lead exposure (up to 10-20 μg/dl) in which the level of brain lead increases linearly, and relatively mildly with the exposure level. If exposure is increased past this level, the brain calcium level exhibits a nonlinear decay, as observed in data (see Fig. 17b, blue curve). If the brain sensitivity to lead is increased (compromised adults or children), the system transitions into a high sensitivity regime (bistability), in which small changes to the intake will have more significant consequences, with very different term outcomes that depend on the history of the system (see Fig. 17b, red and green curves). In this case, the bistability window can be proposed as the mechanism underlying the nonlinear behavior observed in data (with the dependence of behavior on the long-term history of the subjects having been confirmed empirically). Even with a fixed BBB sensitivity, a history of exposure will prime the subject towards the steady state with higher brain lead, and lower calcium (poorer cognitive performance). The transition window occurs in children at much lower levels of lead intake and blood concentrations: since their BBB efficiency and sensitivity are lower, the bistability window shifts to values of p extremely close to zero. This explains why, due to their systemic sensitivity to low level contamination during the developmental stages, children remain a population highly susceptible to lead poisoning, even after great effort has been invested recently in eliminating most higher level lead contamination.
Limitations and future work. The main limitation of the model is its phenomenological nature. We consider the construction of such a model as a first important step towards further developments. Our conclusions can be interpreted as a proof of principle that a nonlinear model based on the lead/calcium competitive dynamics