BACKGROUND
the quantitative modeling of toxicity pathways is a topic of current interest in predictive pharmacology and toxicology  <cit> .
one of its challenges is to integrate omics data with systems biology models for parametric inference and model checking  <cit> .
in a recent paper, wilmes et al.
<cit>  demonstrated a qualitative integration of transcriptomic , proteomic  and metabolomic  data streams to gain a mechanistic understanding of cyclosporine a  toxicity.
csa is an important molecule for immunosuppressive treatment and is used in many post-graft medical protocols  <cit> .
it is mainly metabolized by cyp3a <dig> and cyp3a <dig> and is a substrate of the p-glycoprotein  efflux transporter  <cit> .
however, at high dose it is nephrotoxic and causes damage to the kidney vasculature, glomerulus and proximal tubule  <cit> .
csa is thought to induce oxidative stress at the mitochondrial level, and co-administration of antioxidants with csa appears to mitigate its nephrotoxic effects  <cit> , yet, the precise mechanisms of its toxicity are still unclear.
the nrf <dig> oxidative response pathway is triggered when oxidative stress is sensed by keap- <dig>  resulting in stabilization and nuclear translocation of nrf <dig>  <cit> .
nrf <dig> binds to the antioxidant response element  inducing the transcription of several genes involved in glutathione synthesis and recycling, antioxidant activity, phase ii metabolism and transport  <cit> .
the nrf <dig> response has been shown to be induced in several tissues in response to chemical or physiological stress.
the kidney and particularly the proximal tubule is especially sensitive to oxidative stress.
many nephrotoxins induce nrf <dig> nuclear translocation and nrf2-dependent gene induction in renal epithelial cells, including potassium bromate, cadmium chloride, diquat dibromide and cyclosporine a  <cit> .
moreover, we have recently shown that physiological stress such as glucose depletion and subsequent re-introduction results in nrf <dig> activation in renal cells  <cit> .
here, we use a subset of the wilmes’ et al.
<cit>  omics data to calibrate the parameters of a systems biology model describing the nrf <dig> pathway.
the model predictive ability is then assessed by comparison to csa toxicity data on rptec cells.
methods
data
rptecs culture conditions, csa concentrations measurements, and tcx, ptx and mtx data collection and analysis were described in detail in wilmes et al.
<cit> .
briefly, rptecs cells were cultured in 3 ml of serum-free medium and matured for two weeks on microporous supports.
they were then treated for fourteen days with daily doses of csa.
the assay medium was renewed prior to each dosing.
three groups of assays were performed in triplicate: control, low csa concentration  and high csa concentration .
csa concentration was measured in the medium on the first day at  <dig>  h, 1 h, 3 h, 6 h and 24 h , on the third, fifth, seventh, and tenth day at 24 h , and on day fourteenth at the same times than on the day one.
intracellular  csa concentration and quantity bound to plastic were measured on the first and last days at the same times.
samples for tcx , ptx  and mtx  measurements were collected at the end of day  <dig>  day  <dig> and day  <dig>  all fold-changes were calculated using the absolute value measured at the first time of the control experiment as a reference, and for all doses.
typical rptec cell volume was determined by electron microscopy and stereology to be 2000 ± 140 μm <dig>  on average,  <dig>  millions cells were present in each assay well.
all those data are given as additional material .
mathematical models
modeling was done into two steps:  modeling the in vitro pharmacokinetics  of csa  with a minimal distribution model.
modeling the effects of csa on omics markers at the cellular level with a coupled pk-systems biology model.
that model was calibrated by markov chain monte carlo  sampling  <cit>  with the above in vitro pk and omics data used together.
calibration summarizes and integrates the information brought by the various types of omics data into the coherent scheme of a unified model.
the calibrated model was then run to predict various quantities of interest at a higher level in the hierarchy of biological effects.
such predictions can be compared to further observations, for example on toxicity, to help check the model.
that process is shown schematically on figure  <dig>
in vitro pharmacokinetic model.
a 3-compartment model was developed to describe csa exchange between cell medium, cells and vial walls  <cit> .
in that model, csa can enter and exit the cells, bind to and unbind from the plastic walls and can be metabolized within cells.
several mathematical forms for exchange rates were tested.
the best fit was obtained using a first order entry into cells with michaelis-menten  exit rate, a first order attachment to vial wall with non-integer  order detachment, and michaelis-menten metabolism.
the following differential equations were used to describe the time course of csa quantities in the cytosol, medium, and on vial walls:
∂csacytosol∂t=clin1csaextracellularvextracellular−clout1⋅csacytosolvcytosol⋅kmout1+csacytosol−vmax⋅csacytosolkm2+csacytosol
∂csaextracellular∂t=−clin1csaextracellularvextracellular+clout1⋅csacytosolvcytosol⋅kmout1+csacytosol−k1⋅csaextracellular+k2csawallk <dig>
∂csawall∂t=k1⋅csaextracellular−k2csawallk <dig>
the model parameters are described in table  <dig>
in vitro
csa kinetic parameters description and their statistical distributions
v
max
*: lu: log-uniform distribution .
**:  <dig> zmol =  <dig> zeptomole = 10−21 mol.
coupled pk-systems biology model of the nrf <dig> pathways .
the model used was adapted from zhang et al.
<cit> .
the model full set of equations is provided in the additional file 1: section  <dig>  in brief, csa induces oxidative stress by increasing reactive oxygen species  production.
ros, owing to their electrophilicity, can be detected by the molecular sensor kelch-like ech-associated protein  <dig> , which promotes the ubiquitination and eventual degradation of nrf <dig>  <cit> .
when keap <dig> is oxidized, nrf <dig> ubiquitination is lowered  <cit> , making nrf <dig> available to enter the nucleus.
once in the nucleus, nrf <dig> binds to small maf proteins to form nrf2-maf heterodimers  <cit> .
those can bind to antioxidant responses elements  in the promoter region of glutamate cysteine ligase catalytic subunit , glutamate cysteine ligase modifier subunit , glutathione synthetase , glutathione peroxidase , and abcc <dig> genes  <cit> , maher,  <dig> # <dig>  gclc, gclm, and gs are involved in gsh synthesis.
gpx detoxifies ros, using gsh as a co-substrate.
zhang’s model was developed for a generic xenobiotic, so the following structural changes were made to consistently describe the cell kinetics and mode of action of csa:
– csa can enter or exit the cell, and attach to or detach from the vial walls as in the in vitro pk model .
inside the cell, csa distribution to the nucleus is also modeled ;
– in the cell, csa is metabolized by cytochrome p <dig> 3a <dig>  into a metabolite x’ .
csa is mainly metabolized by cyp3a isoforms  <cit> , and in kidney cells in vivo only cyp3a <dig> is significantly expressed  <cit> ;
– oxidative stress, characterized by the total quantity of oxidative compounds in the cell , was explicitly introduced in the model as a state variable ;
– the production of ros depends on csa concentration in the cell, and ros are eliminated by gpx  in a non-reactive species pool  ;
– keap <dig> and the nrf2-keap <dig> complex are oxidized by ros .
all the other equations are the same as in zhang et al.
<cit> .
in addition, some parameters had to be set to particular values for csa: in the model, xenobiotics can bind to the ahr nuclear receptor.
however, csa is not a known ahr ligand, so its binding parameters  were set to zero.
additional file 1: table s <dig> gives the model parameters’ values and the state variables’ initial values.
calibration of the models
bayesian inference  <cit>  was used to calibrate all the in vitro pk model parameters with data on the csa quantities measured in the medium, cells, and on vial walls .
non-informative  prior parameter distributions were used .
their bounds were set to span an arbitrarily large range of values without constraining estimation.
the data likelihoods were assumed to follow a lognormal distribution around the model predictions, a standard assumption with such measurements.
measurement errors’ geometric standard deviations were assumed to be specific of each of the three measurement types .
they were assigned vague log-uniform prior distributions  and were calibrated together with the other model parameters.
therefore a total of  <dig> parameters  were calibrated.
the posterior statistical distributions of those parameters were obtained by mcmc sampling  <cit> .
for the coupled pk-nrf <dig> pathway model, the parameters directly controlling csa kinetics were set to the joint posterior distribution mode found by the above calibration .
another  <dig> structural parameters  were selected for calibration on the basis of our understanding of the model structure with the help of a preliminary monte carlo sensitivity analysis  <cit>  .
they were calibrated using fold-change omics data  on nrf <dig> mrna, cyp3a <dig> mrna, gs mrna, gclc mrna, gclm mrna, gst mrna, gpx mrna, abcc <dig> mrna, gclm protein, gs protein, mrp <dig> protein, γ-glutamylcysteine , and gsh.
four of the parameters calibrated have a direct influence on the rate of ros synthesis, metabolism and interaction with keap <dig>  another  <dig> parameters control the activation and induction of nrf <dig>  gclc, gclm, gst, gpx, cyp3a <dig>  gs and abcc <dig> genes transcription.
another six parameters control synthesis and degradation of γ-gc and gsh, and the two last parameters control csa metabolism and nrf <dig> and maf binding.
model predicted fold-changes were computed the same way as the experimental ones, using the quantity predicted at the first time of the control experiment as a reference.
the prior parameter distributions chosen were either vague uniform distributions  or lognormal distributions centered around the values used by zhang  <cit>  with a geometric sd corresponding to a factor  <dig> .
the data likelihoods were assumed to be lognormal distributions around the model predictions.
the same measurement error geometric standard deviation was assumed for all omics measurements.
it was calibrated together with the other model parameters, using a vague log-uniform prior.
here also, posterior distributions of the parameter values were obtained by mcmc sampling.
for each model parameter sampled, convergence was evaluated by computing the potential scale reduction criterion of gelman and rubin  <cit> .
k
ros
*:  <dig> zmol =  <dig> zeptomole = 10−21 mol.
**: ln: log-normal distribution .
quantification of csa toxicity for rptecs
the csa concentration causing a 50% decrease  in rptecs’ viability was estimated from the data of limonciel et al.
<cit>  who report dose–response data on the viability of rptecs, 3 t <dig> and heparg cells after exposure to various chemicals, including csa.
the dose range for csa was not large enough to directly estimate an ec <dig>  but a dose–response relationship could nevertheless be reconstructed by meta-analysis in a bayesian framework, borrowing information from the full dose–response observed in the more sensitive 3 t <dig> and heparg cells.
a standard decreasing log-logistic model  <cit>  was calibrated to those dose-viability data using mcmc simulations .
software used
all model simulations and mcmc calibrations were performed with gnu mcsim v <dig> .
<dig>  <cit> .
the r software, version  <dig> .
<dig>  <cit>  was used for other statistical analyses and plots.
RESULTS
in vitro pharmacokinetic model
results for the in vitro pharmacokinetic model have been previously reported in wilmes et al.
<cit>  and are briefly summarized here.
overall, the data were well simulated .
exposure to low concentrations of csa  led to a dynamic steady state in about 3 days.
the average ratio between extracellular and intracellular csa concentrations was about  <dig>  consistent with the fact that csa is lipophilic and accumulates in cells.
exposure to high concentrations of csa  led to major alterations of the biodistribution of csa over time.
steady state was reached in the cells only after approximately 7 days.
the average ratio of intracellular to extracellular csa concentrations was about  <dig> on the first day  and around  <dig> on the last day.
those results highlight the presence of nonlinear phenomena in the distribution of csa in rptecs.
note here the importance of administering repeated doses: a unique administration would not have uncovered those phenomena.
figure  <dig> show the influence of the extra-cellular concentration of csa on the evolution of intra-cellular csa quantities over time.
above about 15 μm extra-cellular csa, intra-cellular concentration does not reach a plateau within 14 days.
coupled pk-systems biology model of the nrf <dig> pathway
all the analyses and predictions presented here were made using a  random sample of  <dig> parameter vectors, obtained by keeping one in each  <dig> of the second half of  <dig>  iterations of five mcmc chains.
approximate convergence required that number of runs, and the median of the gelman and rubin’s criterion was  <dig> .
figures  <dig> and  <dig> show the model fit obtained for the in vitro omics data, at low and high csa exposure dose, respectively.
the bundle of curves presented was obtained using the maximum posterior parameter vector and  <dig> other parameter sets randomly drawn from their joint posterior distribution.
it reflects residual uncertainty in the model predictions, resulting from unavoidable measurement errors and modeling approximations.
overall the data  are quite well fitted.
after calibration, the relative differences between data and predictions are 38% on average .
that is somewhat above, but not by much, the expected measurement precision of the data.
the worst fits are seen for the genes whose transcription apparently decreased, while the model could only predict an increased transcription.
yet, in most cases the decrease observed was modest .
also, the early and rather large increase in glutathione peroxidase mrna at 5 μm csa, with a decrease at 15 μm, could not be reproduced.
for all species, the time profiles are clearly different between the two doses.
at 5 μm csa exposure  periodic oscillations are pervasive.
for many curves  the system does not appear to have reached a dynamic equilibrium within 14 days.
the oscillations’ period may differ from one curve to another and goes up to four days, even though the period of csa administration was exactly one day.
additional file 1: figure s <dig> extends the simulation length to 60 days, time at which a dynamic equilibrium is reached in all cases with the same oscillation pattern.
at 15 μm csa exposure  two patterns emerge: the first type of profile, which concerns all species except gsh and γ-gc, is a plateauing curve.
different maximum values are reached after three days by different curves.
the second type, for gsh and γ-gc, starts with oscillations which do not stabilized within 14 days.
additional file 1: figure s <dig> extends the simulation length to 60 days and shows more clearly the long-term behavior of gsh and γ-gc.
the initial oscillations decrease gradually in amplitude and completely disappear after about 30 days.
figure  <dig> shows model predictions of the time course for the quantity of two non-observed chemical species – the free nuclear nrf <dig> protein and cellular ros – over 14 days, at either low or high repeated csa dosing.
additional simulations were performed up to 60 days and the trends were similar .
as for the previous chemical species for which we had data, large differences are seen between low and high dosing.
at low csa exposure a cyclic pattern is observed, which disappears at high exposure where the ros quantity grows  while the nrf <dig> protein quantity systematically reaches a plateau.
figure  <dig> shows 3d plots of the influence of the extracellular csa concentration on the time course of cellular ros, nuclear nrf <dig> protein, cellular gsh and cellular gcl quantities.
figures for other species  are not shown because their profiles are very similar to the gcl one.
the extracellular concentration of csa has a large influence on the amount of ros in the cytosol.
for extracellular csa concentrations below 8 μm csa, the concentration  vs. time profile of cytosolic ros is oscillating, above 8 μm csa, the ros profile rises in a hockey-stick fashion.
for nuclear nrf <dig>  cellular gsh and cellular gcl, depending on the extracellular csa concentration, the model predicts either oscillating or plateauing profiles.
as for ros, the transition is rather abrupt and occurs approximately at 8 μm csa.
comparison to csa ec <dig> in rptecs
based on our laboratory data, an estimate of  <dig>  μm  was obtained for the ec <dig> of csa for its effect on rptecs viability , after 14 days of repeated dosing.
note that the maximum concentration for csa in water is about 50 μm so that value would actually be difficult to exceed.
in any case, the above range of csa ec <dig> in rptecs matches the predicted increase in ros beyond 15 μm seen on figure  <dig>
discussion
a proper assessment of drug or chemical safety from in vitro assays requires the measurement of concentration of the parent molecule  in the assay medium and in cells  <cit> .
joint kinetic and effect modeling can then be used to interpolate and extrapolate the data obtained.
here, an in vitro pharmacokinetic model was first built using lc-ms/ms data on the distribution of csa over time in human rptecs.
it was then extended to include a description of the nrf <dig> pathway response to the resulting oxidative stress.
csa is highly lipophilic and its rapid uptake and accumulation in cells was observed.
at 5 μm csa , the model indicated that steady-state was reached in about 2 days, whereas at 15 μm csa, steady-state was reached only after 7 days.
moreover, cellular csa concentrations at steady-state were clearly not proportional to exposure, and a disproportionate accumulation of csa was observed at high exposure.
that could be explained by an interplay between the saturation of csa metabolism and transport by p-glycoprotein out of the cells.
however, the amount of csa metabolized in our in vitro pk experiments on rptecs seems limited to only about 15% of the total dose applied, at either low or high dose.
so metabolism alone cannot explain the large increase in concentration that was observed.
also, if csa entry and exit from rptecs were simply linear, the ratio of intra-cellular to extra-cellular concentrations would stay constant in time, and both curves would be parallel on the log-scale .
however, the concentration ratio clearly increases with time and that is easily explained by the saturation of the p-glycoprotein efflux mechanism which was observed by wilmes et al.
<cit>  above 5 μm csa exposure.
drug accumulation in target tissues is often associated with tissue-specific toxicities, and it is important to account for it.
however, we did not observe a direct modulation of csa pk by its pd in our in vitro system, even though csa interactions with transporters are known  <cit> .
in particular, cyp 3a <dig> levels were not affected by csa levels, so csa metabolism was not disturbed by induction or repression.
zhang’s nrf <dig> model was not intended to be used specifically with csa or our cell system, so we had to re-calibrate several parameter values.
this was done in a bayesian statistical framework  <cit> , to take into account the prior information we had on several parameters.
convergence of the mcmc simulations was difficult to obtain due to the high nonlinearity of the model and the presence of cycles.
basically almost any sub-multiple of the actual period would gives an acceptable fit, given the measurement noise.
fortunately, the prior parameter values documented by zhang et al.
stabilize the estimation process.
it would be probably impossible to calibrate the model with a simple maximum likelihood approach .
for most parameters, the posterior mean estimate was clearly different from the prior mean.
the parameters controlling ros formation  and ros elimination  had respectively higher and lower values, compared to zhang’s model, after mcmc sampling.
we centered our gpx parameters’ priors on the values used by zhang et al.
for gst.
since the observed gst and gpx transcriptomic data profiles were very different, it is not surprising that the posterior distributions of gpx parameters are clearly different from their prior.
via the nrf <dig> pathway, csa seems to have a large influence on glutamate cysteine ligase synthesis.
while the basal transcription rates of gclc and gclm have posterior values close to those of zhang’s et al., parameters of gclc and gclm genes regulation by nrf <dig>  linked to ros and csa levels, have at values about four times higher.
indeed the model is imperfect in that it does not describe many additional controls and cross-talks with other pathways, or makes approximations.
for example, gpx is specific for h2o <dig>  and ros are eliminated by a number of enzymes, some of which nrf <dig> influences negatively.
consequently, for example, our model cannot explain the  decreases that were observed for some gene transcripts and does not describe well the time course of gpx mrna at low dose.
more observations or model refinements would be needed to understand the origin  of that discrepancy.
the model still gives access to unmeasured effects of csa to cells, closer to a toxicity endpoint.
the generation of ros by csa is an important toxicity mechanism for that molecule.
the retro-control of ros scavenging by the nrf <dig> signaling pathway induces a highly nonlinear behavior illustrated on figure  <dig>  the cyclic patterns observed at low csa exposure  are interesting.
a recent paper describes circadian oscillations of the nrf2/gsh pathway in mice lung  <cit> .
that pathway is well conserved and present in most cells since it regulates oxidative species generated during respiration.
therefore, a rhythmic pattern of nrf <dig> activity in rptecs would not be surprising even in the absence of csa.
in addition, the daily exposure of rptecs to 5 μm csa is likely to have produced a manageable burst of oxidative stress at the beginning of each day.
that alone could explain the cycles seen, even if the nonlinear dynamics of the model result in a period of two to three days for those cycles.
ros generation runs out of control at csa exposure levels close to the high dose assayed in vitro  and the cycles disappear .
we have an external corroboration of this finding: the 15 μm concentration was experimentally chosen to be the highest not affecting cell survival.
above that level, toxicity starts to have an impact on survival and the in vitro ec <dig> of csa in rptecs is probably close to 40 μm, so our model predictions seem reasonable.
however, as in many systems biology models, only one signaling pathway has been taken into account.
we also do not have extensive data allowing for an in-depth statistical cross-validation of the many components of the model.
other ros scavenging mechanisms are present in rptecs and could be involved during csa exposure.
on the other hand, csa nephrotoxicity involves several mechanism  <cit>  and it is possible that ros generation is not alone in causing critical damages.
CONCLUSIONS
integrating omics approaches with mathematical systems biology models is still rarely done  <cit> , even though that seems the best way to both understand the data and improve the predictive ability of the models  <cit> .
our modeling and simulations of the csa mediated ros production gives biologists insight into mechanisms of toxicity and provide quantitative estimates of toxicity beyond the time and dose range used in experiments.
to go further, it would be interesting to have a more precise model description of gsh synthesis in the model, since cellular ros concentrations are clearly correlated to gsh.
it would also be interesting to couple this model with a physiological based pharmacokinetic  model for csa to be able to better predict human response.
still, our results demonstrate the possibility to use different omic data streams to extrapolate in time and dose the response of the nrf <dig> pathway to oxidative damage, far beyond our current experimental possibilities.
abbreviations
abcc2: atp-binding cassette sub-family c member 2; ahr: aryl hydrocarbon receptor ; are: antioxidant response element; arnt: aryl hydrocarbon receptor nuclear translocator ; csa: cyclosporine a; csacytosol: cystosol csa quantity; csaextracellular: extracellular csa quantity; csawall: csa on wall quantity; cyp: cytochrome p <dig> 3a5; dre: dioxin response element; gcl: glutamate cysteine ligase; gclc: glutamate cysteine ligase catalytic subunit; gclm: glutamate cysteine ligase modifier subunit; gpx: glutathione peroxidase; gs: glutathione synthetase; gsh: glutathione; gssg: glutathione disulfide; keap1: kelch-like ech-associated protein 1; mcmc: markov chain monte carlo; mrna: messenger ribonucleic acid; mrp2: multidrug resistance associated protein 2; mtx: metabolomic; nma: nrf2-maf-are complex; nrf2: nuclear factor -like 2; nrs: non reactive species; pk: pharmacokinetics; ptx: proteomic; ros: reactive oxygen species; rptec: renal proximal tubular epithelial cells; tcx: transcriptomic; xaa: complex of csa-ahr-arnt; xaad: complex of csa-ahr-arnt-dre; γ-gc: γ-glutamylcysteine.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
jh and fyb coded and ran the model.
pj performed the experiments.
all authors participated in the redaction of this article.
all authors read and approved the final manuscript.
supplementary material
additional file 1
section  <dig>  differential equations of the nrf <dig> pathway model’.
section  <dig>  preliminary sensitivity analysis for the selection of nrf <dig> model parameters to calibrate.
section  <dig>  quantification of csa toxicity for rptecs.
figure s <dig>  maximum posterior fits of the log-logistic viability model for 3t <dig> and heparg cells viability data as a function of csa exposure concentration.
figure s <dig>  maximum posterior fit and 95% confidence bounds of the log-logistic viability model for rptecs viability data as a function of csa exposure concentration.
table s <dig>  model parameters values and initial state variables values.
table s <dig>  cyclosporine a quantities measured in the extracellular medium at low csa concentration exposure.
table s <dig>  intracellular cyclosporine a quantities measured at low csa concentration exposure.
table s <dig>  cyclosporine a quantities measured on plastic at low csa concentration exposure.
table s <dig>  cyclosporine a quantities measured in the extracellular medium at high csa concentration exposure.
table s <dig>  intracellular cyclosporine a quantities measured at high csa concentration exposure.
table s <dig>  cyclosporine a quantities measured on plastic at high csa concentration exposure.
table s <dig>  fold changes measured at low csa concentration.
table s <dig>  fold changes measured at high csa concentration.
figure s <dig>  model fit to the data.
the data values are plotted against the model predictions, after model calibration.
the pk data are represented by black circles, the metabolomic data by green square, transcriptomic by red triangles and proteomics by blue inverted triangles.
figure s <dig>  transcriptomics, proteomics, and metabolomics  fold-changes time-course in rptec cells during  <dig> days with repeated low dose csa dosing.
figure s <dig>  transcriptomics, proteomics, and metabolomics  fold-changes time-course in rptec cells during  <dig> days with repeated high dose csa dosing.
click here for file
acknowledgements
this work was supported by the european commission, 7th fp project predict-iv , the nc3r crack-it award , and the french ministry for the environment .

 the quantitative modeling of toxicity pathways is a topic of current interest in predictive pharmacology and toxicology. one of its challenges is to integrate omics data with systems biology models for parametric inference and model checking. in a recent paper, wilmes et al. demonstrated a qualitative integration of transcriptomic, proteomic and metabolomic data streams to gain a mechanistic understanding of cyclosporine a toxicity. csa is an important molecule for immunosuppressive treatment and is used in many post-graft medical protocols. it is mainly metabolized by cyp3a and cyp3a and is a substrate of the p-glycoprotein efflux transporter. however, at high dose it is nephrotoxic and causes damage to the kidney vasculature, glomerulus and proximal tubule. csa is thought to induce oxidative stress at the mitochondrial level, and co-administration of antioxidants with csa appears to mitigate its nephrotoxic effects, yet, the precise mechanisms of its toxicity are still unclear. the nrf oxidative response pathway is triggered when oxidative stress is sensed by keap- resulting in stabilization and nuclear translocation of nrf. nrf binds to the antioxidant response element inducing the transcription of several genes involved in glutathione synthesis and recycling, antioxidant activity, phase ii metabolism and transport. the nrf response has been shown to be induced in several tissues in response to chemical or physiological stress. the kidney and particularly the proximal tubule is especially sensitive to oxidative stress. many nephrotoxins induce nrf nuclear translocation and nrf2-dependent gene induction in renal epithelial cells, including potassium bromate, cadmium chloride, diquat dibromide and cyclosporine a. moreover, we have recently shown that physiological stress such as glucose depletion and subsequent re-introduction results in nrf activation in renal cells. here, we use a subset of the wilmes’ et al. omics data to calibrate the parameters of a systems biology model describing the nrf pathway. the model predictive ability is then assessed by comparison to csa toxicity data on rptec cells. methods data rptecs culture conditions, csa concentrations measurements, and tcx, ptx and mtx data collection and analysis were described in detail in wilmes et al.. briefly, rptecs cells were cultured in 3 ml of serum-free medium and matured for two weeks on microporous supports. they were then treated for fourteen days with daily doses of csa. the assay medium was renewed prior to each dosing. three groups of assays were performed in triplicate: control, low csa concentration and high csa concentration. csa concentration was measured in the medium on the first day at  h, 1 h, 3 h, 6 h and 24 h, on the third, fifth, seventh, and tenth day at 24 h, and on day fourteenth at the same times than on the day one. intracellular csa concentration and quantity bound to plastic were measured on the first and last days at the same times. samples for tcx, ptx and mtx measurements were collected at the end of day day and day all fold-changes were calculated using the absolute value measured at the first time of the control experiment as a reference, and for all doses. typical rptec cell volume was determined by electron microscopy and stereology to be 2000 ± 140 μm on average, millions cells were present in each assay well. all those data are given as additional material. mathematical models modeling was done into two steps: modeling the in vitro pharmacokinetics of csa with a minimal distribution model. modeling the effects of csa on omics markers at the cellular level with a coupled pk-systems biology model. that model was calibrated by markov chain monte carlo sampling with the above in vitro pk and omics data used together. calibration summarizes and integrates the information brought by the various types of omics data into the coherent scheme of a unified model. the calibrated model was then run to predict various quantities of interest at a higher level in the hierarchy of biological effects. such predictions can be compared to further observations, for example on toxicity, to help check the model. that process is shown schematically on figure  in vitro pharmacokinetic model. a 3-compartment model was developed to describe csa exchange between cell medium, cells and vial walls. in that model, csa can enter and exit the cells, bind to and unbind from the plastic walls and can be metabolized within cells. several mathematical forms for exchange rates were tested. the best fit was obtained using a first order entry into cells with michaelis-menten exit rate, a first order attachment to vial wall with non-integer order detachment, and michaelis-menten metabolism. the following differential equations were used to describe the time course of csa quantities in the cytosol, medium, and on vial walls: ∂csacytosol∂t=clin1csaextracellularvextracellular−clout1⋅csacytosolvcytosol⋅kmout1+csacytosol−vmax⋅csacytosolkm2+csacytosol ∂csaextracellular∂t=−clin1csaextracellularvextracellular+clout1⋅csacytosolvcytosol⋅kmout1+csacytosol−k1⋅csaextracellular+k2csawallk ∂csawall∂t=k1⋅csaextracellular−k2csawallk the model parameters are described in table  in vitro csa kinetic parameters description and their statistical distributions v max *: lu: log-uniform distribution. **: zmol =  zeptomole = 10−21 mol. coupled pk-systems biology model of the nrf pathways. the model used was adapted from zhang et al.. the model full set of equations is provided in the additional file 1: section in brief, csa induces oxidative stress by increasing reactive oxygen species production. ros, owing to their electrophilicity, can be detected by the molecular sensor kelch-like ech-associated protein, which promotes the ubiquitination and eventual degradation of nrf. when keap is oxidized, nrf ubiquitination is lowered, making nrf available to enter the nucleus. once in the nucleus, nrf binds to small maf proteins to form nrf2-maf heterodimers. those can bind to antioxidant responses elements in the promoter region of glutamate cysteine ligase catalytic subunit, glutamate cysteine ligase modifier subunit, glutathione synthetase, glutathione peroxidase, and abcc genes, maher, # gclc, gclm, and gs are involved in gsh synthesis. gpx detoxifies ros, using gsh as a co-substrate. zhang’s model was developed for a generic xenobiotic, so the following structural changes were made to consistently describe the cell kinetics and mode of action of csa: – csa can enter or exit the cell, and attach to or detach from the vial walls as in the in vitro pk model. inside the cell, csa distribution to the nucleus is also modeled; – in the cell, csa is metabolized by cytochrome p 3a into a metabolite x’. csa is mainly metabolized by cyp3a isoforms, and in kidney cells in vivo only cyp3a is significantly expressed; – oxidative stress, characterized by the total quantity of oxidative compounds in the cell, was explicitly introduced in the model as a state variable; – the production of ros depends on csa concentration in the cell, and ros are eliminated by gpx in a non-reactive species pool; – keap and the nrf2-keap complex are oxidized by ros. all the other equations are the same as in zhang et al.. in addition, some parameters had to be set to particular values for csa: in the model, xenobiotics can bind to the ahr nuclear receptor. however, csa is not a known ahr ligand, so its binding parameters were set to zero. additional file 1: table s gives the model parameters’ values and the state variables’ initial values. calibration of the models bayesian inference was used to calibrate all the in vitro pk model parameters with data on the csa quantities measured in the medium, cells, and on vial walls. non-informative prior parameter distributions were used. their bounds were set to span an arbitrarily large range of values without constraining estimation. the data likelihoods were assumed to follow a lognormal distribution around the model predictions, a standard assumption with such measurements. measurement errors’ geometric standard deviations were assumed to be specific of each of the three measurement types. they were assigned vague log-uniform prior distributions and were calibrated together with the other model parameters. therefore a total of parameters were calibrated. the posterior statistical distributions of those parameters were obtained by mcmc sampling. for the coupled pk-nrf pathway model, the parameters directly controlling csa kinetics were set to the joint posterior distribution mode found by the above calibration. another structural parameters were selected for calibration on the basis of our understanding of the model structure with the help of a preliminary monte carlo sensitivity analysis. they were calibrated using fold-change omics data on nrf mrna, cyp3a mrna, gs mrna, gclc mrna, gclm mrna, gst mrna, gpx mrna, abcc mrna, gclm protein, gs protein, mrp protein, γ-glutamylcysteine, and gsh. four of the parameters calibrated have a direct influence on the rate of ros synthesis, metabolism and interaction with keap another parameters control the activation and induction of nrf gclc, gclm, gst, gpx, cyp3a gs and abcc genes transcription. another six parameters control synthesis and degradation of γ-gc and gsh, and the two last parameters control csa metabolism and nrf and maf binding. model predicted fold-changes were computed the same way as the experimental ones, using the quantity predicted at the first time of the control experiment as a reference. the prior parameter distributions chosen were either vague uniform distributions or lognormal distributions centered around the values used by zhang with a geometric sd corresponding to a factor. the data likelihoods were assumed to be lognormal distributions around the model predictions. the same measurement error geometric standard deviation was assumed for all omics measurements. it was calibrated together with the other model parameters, using a vague log-uniform prior. here also, posterior distributions of the parameter values were obtained by mcmc sampling. for each model parameter sampled, convergence was evaluated by computing the potential scale reduction criterion of gelman and rubin. k ros *: zmol =  zeptomole = 10−21 mol. **: ln: log-normal distribution. quantification of csa toxicity for rptecs the csa concentration causing a 50% decrease in rptecs’ viability was estimated from the data of limonciel et al. who report dose–response data on the viability of rptecs, 3 t and heparg cells after exposure to various chemicals, including csa. the dose range for csa was not large enough to directly estimate an ec but a dose–response relationship could nevertheless be reconstructed by meta-analysis in a bayesian framework, borrowing information from the full dose–response observed in the more sensitive 3 t and heparg cells. a standard decreasing log-logistic model was calibrated to those dose-viability data using mcmc simulations. software used all model simulations and mcmc calibrations were performed with gnu mcsim v.. the r software, version. was used for other statistical analyses and plots.  in vitro pharmacokinetic model results for the in vitro pharmacokinetic model have been previously reported in wilmes et al. and are briefly summarized here. overall, the data were well simulated. exposure to low concentrations of csa led to a dynamic steady state in about 3 days. the average ratio between extracellular and intracellular csa concentrations was about consistent with the fact that csa is lipophilic and accumulates in cells. exposure to high concentrations of csa led to major alterations of the biodistribution of csa over time. steady state was reached in the cells only after approximately 7 days. the average ratio of intracellular to extracellular csa concentrations was about on the first day and around on the last day. those results highlight the presence of nonlinear phenomena in the distribution of csa in rptecs. note here the importance of administering repeated doses: a unique administration would not have uncovered those phenomena. figure  show the influence of the extra-cellular concentration of csa on the evolution of intra-cellular csa quantities over time. above about 15 μm extra-cellular csa, intra-cellular concentration does not reach a plateau within 14 days. coupled pk-systems biology model of the nrf pathway all the analyses and predictions presented here were made using a random sample of parameter vectors, obtained by keeping one in each of the second half of iterations of five mcmc chains. approximate convergence required that number of runs, and the median of the gelman and rubin’s criterion was. figures  and show the model fit obtained for the in vitro omics data, at low and high csa exposure dose, respectively. the bundle of curves presented was obtained using the maximum posterior parameter vector and other parameter sets randomly drawn from their joint posterior distribution. it reflects residual uncertainty in the model predictions, resulting from unavoidable measurement errors and modeling approximations. overall the data are quite well fitted. after calibration, the relative differences between data and predictions are 38% on average. that is somewhat above, but not by much, the expected measurement precision of the data. the worst fits are seen for the genes whose transcription apparently decreased, while the model could only predict an increased transcription. yet, in most cases the decrease observed was modest. also, the early and rather large increase in glutathione peroxidase mrna at 5 μm csa, with a decrease at 15 μm, could not be reproduced. for all species, the time profiles are clearly different between the two doses. at 5 μm csa exposure periodic oscillations are pervasive. for many curves the system does not appear to have reached a dynamic equilibrium within 14 days. the oscillations’ period may differ from one curve to another and goes up to four days, even though the period of csa administration was exactly one day. additional file 1: figure s extends the simulation length to 60 days, time at which a dynamic equilibrium is reached in all cases with the same oscillation pattern. at 15 μm csa exposure two patterns emerge: the first type of profile, which concerns all species except gsh and γ-gc, is a plateauing curve. different maximum values are reached after three days by different curves. the second type, for gsh and γ-gc, starts with oscillations which do not stabilized within 14 days. additional file 1: figure s extends the simulation length to 60 days and shows more clearly the long-term behavior of gsh and γ-gc. the initial oscillations decrease gradually in amplitude and completely disappear after about 30 days. figure  shows model predictions of the time course for the quantity of two non-observed chemical species – the free nuclear nrf protein and cellular ros – over 14 days, at either low or high repeated csa dosing. additional simulations were performed up to 60 days and the trends were similar. as for the previous chemical species for which we had data, large differences are seen between low and high dosing. at low csa exposure a cyclic pattern is observed, which disappears at high exposure where the ros quantity grows while the nrf protein quantity systematically reaches a plateau. figure  shows 3d plots of the influence of the extracellular csa concentration on the time course of cellular ros, nuclear nrf protein, cellular gsh and cellular gcl quantities. figures for other species are not shown because their profiles are very similar to the gcl one. the extracellular concentration of csa has a large influence on the amount of ros in the cytosol. for extracellular csa concentrations below 8 μm csa, the concentration vs. time profile of cytosolic ros is oscillating, above 8 μm csa, the ros profile rises in a hockey-stick fashion. for nuclear nrf cellular gsh and cellular gcl, depending on the extracellular csa concentration, the model predicts either oscillating or plateauing profiles. as for ros, the transition is rather abrupt and occurs approximately at 8 μm csa. comparison to csa ec in rptecs based on our laboratory data, an estimate of  μm was obtained for the ec of csa for its effect on rptecs viability, after 14 days of repeated dosing. note that the maximum concentration for csa in water is about 50 μm so that value would actually be difficult to exceed. in any case, the above range of csa ec in rptecs matches the predicted increase in ros beyond 15 μm seen on figure  discussion a proper assessment of drug or chemical safety from in vitro assays requires the measurement of concentration of the parent molecule in the assay medium and in cells. joint kinetic and effect modeling can then be used to interpolate and extrapolate the data obtained. here, an in vitro pharmacokinetic model was first built using lc-ms/ms data on the distribution of csa over time in human rptecs. it was then extended to include a description of the nrf pathway response to the resulting oxidative stress. csa is highly lipophilic and its rapid uptake and accumulation in cells was observed. at 5 μm csa, the model indicated that steady-state was reached in about 2 days, whereas at 15 μm csa, steady-state was reached only after 7 days. moreover, cellular csa concentrations at steady-state were clearly not proportional to exposure, and a disproportionate accumulation of csa was observed at high exposure. that could be explained by an interplay between the saturation of csa metabolism and transport by p-glycoprotein out of the cells. however, the amount of csa metabolized in our in vitro pk experiments on rptecs seems limited to only about 15% of the total dose applied, at either low or high dose. so metabolism alone cannot explain the large increase in concentration that was observed. also, if csa entry and exit from rptecs were simply linear, the ratio of intra-cellular to extra-cellular concentrations would stay constant in time, and both curves would be parallel on the log-scale. however, the concentration ratio clearly increases with time and that is easily explained by the saturation of the p-glycoprotein efflux mechanism which was observed by wilmes et al. above 5 μm csa exposure. drug accumulation in target tissues is often associated with tissue-specific toxicities, and it is important to account for it. however, we did not observe a direct modulation of csa pk by its pd in our in vitro system, even though csa interactions with transporters are known. in particular, cyp 3a levels were not affected by csa levels, so csa metabolism was not disturbed by induction or repression. zhang’s nrf model was not intended to be used specifically with csa or our cell system, so we had to re-calibrate several parameter values. this was done in a bayesian statistical framework, to take into account the prior information we had on several parameters. convergence of the mcmc simulations was difficult to obtain due to the high nonlinearity of the model and the presence of cycles. basically almost any sub-multiple of the actual period would gives an acceptable fit, given the measurement noise. fortunately, the prior parameter values documented by zhang et al. stabilize the estimation process. it would be probably impossible to calibrate the model with a simple maximum likelihood approach. for most parameters, the posterior mean estimate was clearly different from the prior mean. the parameters controlling ros formation and ros elimination had respectively higher and lower values, compared to zhang’s model, after mcmc sampling. we centered our gpx parameters’ priors on the values used by zhang et al. for gst. since the observed gst and gpx transcriptomic data profiles were very different, it is not surprising that the posterior distributions of gpx parameters are clearly different from their prior. via the nrf pathway, csa seems to have a large influence on glutamate cysteine ligase synthesis. while the basal transcription rates of gclc and gclm have posterior values close to those of zhang’s et al., parameters of gclc and gclm genes regulation by nrf linked to ros and csa levels, have at values about four times higher. indeed the model is imperfect in that it does not describe many additional controls and cross-talks with other pathways, or makes approximations. for example, gpx is specific for h2o and ros are eliminated by a number of enzymes, some of which nrf influences negatively. consequently, for example, our model cannot explain the decreases that were observed for some gene transcripts and does not describe well the time course of gpx mrna at low dose. more observations or model refinements would be needed to understand the origin of that discrepancy. the model still gives access to unmeasured effects of csa to cells, closer to a toxicity endpoint. the generation of ros by csa is an important toxicity mechanism for that molecule. the retro-control of ros scavenging by the nrf signaling pathway induces a highly nonlinear behavior illustrated on figure  the cyclic patterns observed at low csa exposure are interesting. a recent paper describes circadian oscillations of the nrf2/gsh pathway in mice lung. that pathway is well conserved and present in most cells since it regulates oxidative species generated during respiration. therefore, a rhythmic pattern of nrf activity in rptecs would not be surprising even in the absence of csa. in addition, the daily exposure of rptecs to 5 μm csa is likely to have produced a manageable burst of oxidative stress at the beginning of each day. that alone could explain the cycles seen, even if the nonlinear dynamics of the model result in a period of two to three days for those cycles. ros generation runs out of control at csa exposure levels close to the high dose assayed in vitro and the cycles disappear. we have an external corroboration of this finding: the 15 μm concentration was experimentally chosen to be the highest not affecting cell survival. above that level, toxicity starts to have an impact on survival and the in vitro ec of csa in rptecs is probably close to 40 μm, so our model predictions seem reasonable. however, as in many systems biology models, only one signaling pathway has been taken into account. we also do not have extensive data allowing for an in-depth statistical cross-validation of the many components of the model. other ros scavenging mechanisms are present in rptecs and could be involved during csa exposure. on the other hand, csa nephrotoxicity involves several mechanism and it is possible that ros generation is not alone in causing critical damages.  integrating omics approaches with mathematical systems biology models is still rarely done, even though that seems the best way to both understand the data and improve the predictive ability of the models. our modeling and simulations of the csa mediated ros production gives biologists insight into mechanisms of toxicity and provide quantitative estimates of toxicity beyond the time and dose range used in experiments. to go further, it would be interesting to have a more precise model description of gsh synthesis in the model, since cellular ros concentrations are clearly correlated to gsh. it would also be interesting to couple this model with a physiological based pharmacokinetic model for csa to be able to better predict human response. still, our results demonstrate the possibility to use different omic data streams to extrapolate in time and dose the response of the nrf pathway to oxidative damage, far beyond our current experimental possibilities. abbreviations abcc2: atp-binding cassette sub-family c member 2; ahr: aryl hydrocarbon receptor; are: antioxidant response element; arnt: aryl hydrocarbon receptor nuclear translocator; csa: cyclosporine a; csacytosol: cystosol csa quantity; csaextracellular: extracellular csa quantity; csawall: csa on wall quantity; cyp: cytochrome p 3a5; dre: dioxin response element; gcl: glutamate cysteine ligase; gclc: glutamate cysteine ligase catalytic subunit; gclm: glutamate cysteine ligase modifier subunit; gpx: glutathione peroxidase; gs: glutathione synthetase; gsh: glutathione; gssg: glutathione disulfide; keap1: kelch-like ech-associated protein 1; mcmc: markov chain monte carlo; mrna: messenger ribonucleic acid; mrp2: multidrug resistance associated protein 2; mtx: metabolomic; nma: nrf2-maf-are complex; nrf2: nuclear factor -like 2; nrs: non reactive species; pk: pharmacokinetics; ptx: proteomic; ros: reactive oxygen species; rptec: renal proximal tubular epithelial cells; tcx: transcriptomic; xaa: complex of csa-ahr-arnt; xaad: complex of csa-ahr-arnt-dre; γ-gc: γ-glutamylcysteine. competing interests the authors declare that they have no competing interests. authors’ contributions jh and fyb coded and ran the model. pj performed the experiments. all authors participated in the redaction of this article. all authors read and approved the final manuscript. supplementary material additional file 1 section differential equations of the nrf pathway model’. section preliminary sensitivity analysis for the selection of nrf model parameters to calibrate. section quantification of csa toxicity for rptecs. figure s maximum posterior fits of the log-logistic viability model for 3t and heparg cells viability data as a function of csa exposure concentration. figure s maximum posterior fit and 95% confidence bounds of the log-logistic viability model for rptecs viability data as a function of csa exposure concentration. table s model parameters values and initial state variables values. table s cyclosporine a quantities measured in the extracellular medium at low csa concentration exposure. table s intracellular cyclosporine a quantities measured at low csa concentration exposure. table s cyclosporine a quantities measured on plastic at low csa concentration exposure. table s cyclosporine a quantities measured in the extracellular medium at high csa concentration exposure. table s intracellular cyclosporine a quantities measured at high csa concentration exposure. table s cyclosporine a quantities measured on plastic at high csa concentration exposure. table s fold changes measured at low csa concentration. table s fold changes measured at high csa concentration. figure s model fit to the data. the data values are plotted against the model predictions, after model calibration. the pk data are represented by black circles, the metabolomic data by green square, transcriptomic by red triangles and proteomics by blue inverted triangles. figure s transcriptomics, proteomics, and metabolomics fold-changes time-course in rptec cells during days with repeated low dose csa dosing. figure s transcriptomics, proteomics, and metabolomics fold-changes time-course in rptec cells during days with repeated high dose csa dosing. click here for file acknowledgements this work was supported by the european commission, 7th fp project predict-iv, the nc3r crack-it award, and the french ministry for the environment.
posterior statistical distributions of the model parameter values were obtained by markov chain monte carlo sampling in a bayesian framework. we present here a quantitative calibration of a differential equation model of the nrf pathway with a subset of the omics data we collected. omics data can be used to that effect in a bayesian statistical framework which retains prior information about the likely parameter values. the model proposed can be used to analyze and predict cellular response to oxidative stress, provided sufficient data to set its parameters to cell-specific values. bayesian data analysiscyclosporine aglutathione pathwayintegrated omicsnrf pathwayoxidative stressrenal proximal tubule epithelial cellssystems biology in vitro pharmacokinetic data on csa exchange between cells, culture medium and vial walls, and data on the time course of omics markers in response to csa exposure were reasonably well fitted with a coupled pk-systems biology model. incorporation of omic data streams for building improved systems biology models has great potential for improving their predictions of biological outcomes. a complex cyclic pattern of ros production and control emerged at 5 μm csa repeated exposure. we have recently shown that cyclosporine a strongly activates the nuclear factor -like pathway in renal proximal tubular epithelial cells exposed in vitro. shortly above those exposure levels, the model predicts a disproportionate increase in cellular ros quantity which is consistent with an in vitro ec of about 40 μm for csa in rptecs.
 incorporation of omic data streams for building improved systems biology models has great potential for improving their predictions of biological outcomes. we have recently shown that cyclosporine a strongly activates the nuclear factor -like pathway in renal proximal tubular epithelial cells exposed in vitro. we present here a quantitative calibration of a differential equation model of the nrf pathway with a subset of the omics data we collected.  in vitro pharmacokinetic data on csa exchange between cells, culture medium and vial walls, and data on the time course of omics markers in response to csa exposure were reasonably well fitted with a coupled pk-systems biology model. posterior statistical distributions of the model parameter values were obtained by markov chain monte carlo sampling in a bayesian framework. a complex cyclic pattern of ros production and control emerged at 5 μm csa repeated exposure. plateau responses were found at 15 μm exposures. shortly above those exposure levels, the model predicts a disproportionate increase in cellular ros quantity which is consistent with an in vitro ec of about 40 μm for csa in rptecs.  the model proposed can be used to analyze and predict cellular response to oxidative stress, provided sufficient data to set its parameters to cell-specific values. omics data can be used to that effect in a bayesian statistical framework which retains prior information about the likely parameter values. bayesian data analysiscyclosporine aglutathione pathwayintegrated omicsnrf pathwayoxidative stressrenal proximal tubule epithelial cellssystems biology 