Integration of Omics Data and Systems Biology Modeling: Effect of Cyclosporine A on the Nrf2 Pathway in Human Renal Kidneys Cells

In a recent paper, Wilmes et al. demonstrated a qualitative integration of omics data streams to gain a mechanistic understanding of cyclosporine A toxicity. One of their major conclusions was that cyclosporine A strongly activates the nuclear factor (erythroid-derived 2)-like 2 pathway (Nrf2) in renal proximal tubular epithelial cells exposed in vitro. We pursue here the analysis of those data with a quantitative integration of omics data with a differential equation model of the Nrf2 pathway. That was done in two steps: (i) Modeling the in vitro pharmacokinetics of cyclosporine A (exchange between cells, culture medium and vial walls) with a minimal distribution model. (ii) Modeling the time course of omics markers in response to cyclosporine A exposure at the cell level with a coupled PK-systems biology model. Posterior statistical distributions of the parameter values were obtained by Markov chain Monte Carlo sampling. Data were well simulated, and the known in vitro toxic effect EC50 was well matched by model predictions. The integration of in vitro pharmacokinetics and systems biology modeling gives us a quantitative insight into mechanisms of cyclosporine A oxidative-stress induction, and a way to predict such a stress for a variety of exposure conditions.


Introduction
The quantitative modeling of toxicity pathways is a topic of current interest in predictive pharmacology and toxicology [2,3,4]. One of its challenges is to integrate omics data with systemic biology models for parametric inference and model checking [5]. In a recent paper, Wilmes et al.
[1] demonstrated a qualitative integration of transcriptomic (TCX), proteomic (PTX) and metabolomic (MTX) data streams to gain a mechanistic understanding of cyclosporine A (CsA) toxicity. CsA is an important molecule for immunosuppressive treatment and is used in many postgraft medical protocols [6]. However, at high dose it is nephrotoxic, causing damage to the kidney vasculature, glomerulus and proximal tubule [7,8,9]. Yet, the precise mechanisms of its toxicity are still unclear: CsA is thought to induce oxidative stress at the mitochondrial level, and coadministration of antioxidants with CsA appears to mitigate its nephrotoxic effects [10].
The Nrf2 oxidative response pathway is triggered when oxidative stress is sensed by Keap-1, resulting in stabilization and nuclear translocation of Nrf2 [11]. Nrf2 binds to the antioxidant response element (ARE) inducing the transcription of several genes involved in glutathione synthesis and recycling, antioxidant activity and phase II metabolism and transport [11]. The Nrf2 response has been shown to induced in several tissues in response to chemical and physiological stressors. The kidney and particularly the proximal tubule is especially sensitive to oxidative stress and many nephrotoxins induce Nrf2 nuclear translocation and Nrf2 dependent gene induction in renal epithelial cells, including potassium bromate, cadmium chloride, diquat dibromide and cyclosporine A [1,12,13]. Moreover, we have recently shown physiological stress such as glucose depletion and subsequent re-introduction results in Nrf2 activation in renal cells [14]. Here, we pursue here the analysis of that data set with a quantitative integration of omics data with a systems biology model of the Nrf2 pathway. The model predictive ability is then assessed.

Data
RPTECs culture conditions, CsA concentrations measurements, and TCX, PTX and MTX data collection and analysis were described in detail in Wilmes et al. [1]. Briefly, RPTECs cells were cultured in 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 (5 μM) and high CsA concentration (15 μM).
CsA concentration was measured in the medium on the first day at 0.5h, 1h, 3h, 6h and 24h (just before changing the medium), on the third, fifth, seventh, and tenth day at 24h (before changing the medium), and on day fourteenth at the same times than on the day one. Intracellular (cell lysate) CsA concentration and quantity bound to plastic were measured on the first and last days at the same times. Samples for TCX (conducted with Illumina ® HT 12 v3 BeadChip arrays), PTX (conducted with HPLC-MS) and MTX (conducted with direct infusion MS) measurements were obtained at the end of day 1, day 3 and day 14. 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 3 .
All the data used are given in Supplementary tables S1 to S8.

Mathematical model
Modeling was done into two steps: (i) Modeling the in vitro pharmacokinetics (PK) of CsA (exchange between cells, medium and vial walls) with a minimal distribution model. (ii) Modeling the effects of CsA on omics markers as the cellular level with a coupled PK-systems biology model. In vitro pharmacokinetic model. A 3-compartment model was developed to describe CsA exchange between cell medium, cells and vial walls [1]. 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 (saturable) exit rate, a first order attachment to vial wall with non-integer (fractal) 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: The model parameters are described in Table 1.
Coupled PK-systems biology model of the Nrf2 pathways ( Figure 1). The model used was adapted from Zhang et al. [15]. The model full set of equations is provided as supplemental material. In brief, CsA induces oxidative stress by increasing reactive oxygen species (ROS) production. ROS, owing to their electrophilicity, can be detected by the molecular sensor Kelch-like ECH-associated protein 1 (Keap1), which promotes the ubiquitination and eventual degradation of Nrf2 [16,17]. When Keap1 is oxidized, Nrf2 ubiquitination is lowered [17], making Nrf2 available to enter the nucleus. Once in the nucleus, Nrf2 binds to small Maf proteins to form Nrf2-Maf heterodimers [18]. Those can bind to antioxidant responses elements (ARE) in the promoter region of glutamate cysteine ligase catalytic subunit (GCLC), glutamate cysteine ligase modifier subunit (GCLM), glutathione synthetase (GS), glutathione peroxidase (GPx), and MRP genes, inducing their transcriptions [16,18]. 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 (eqs. A10, A11, A13). Inside the cell, CsA distribution to the nucleus is explicitly modeled (eq. A12); -In the cell, CsA is metabolized by cytochrome P450 3A5 (CYP3A5) into a metabolite X′ (not followed and non-influent on the system). CsA is mainly metabolized by CYP3A isoforms [6], and in kidney cells only CYP3A5 is significantly expressed [19]; -Oxidative stress, characterized by the total quantity of oxidative compounds in the cell (ROS), was explicitly introduced in the model as a state variable (eq. A75); -The production of ROS depends on CsA concentration in the cell, and ROS are eliminated by GPx (eq. A75) in a non-reactive species pool (NRS) (not followed and non-influent on the system); -Keap1 and the Nrf2-Keap1 complex are oxidized by ROS (eqs. A52, A53, A72, A73).
All the other equations are the same than in Zhang et al. [15]. 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 (k b 2 and k b 5 ) were set to zero. Supplemental Table S9 gives the model parameters and the state variables initial values.

Calibration of the models
The in vitro PK model parameters were calibrated, using Bayesian inference [5,20], with data on the CsA quantities measured in the medium, cells, and on vial walls (Table S1 to S6). Noninformative (vague) prior parameter distributions were used (see Table 1). The data likelihoods were assumed to follow a lognormal distribution around the model predictions, a standard assumption with such measurements. Measurement error geometric standard deviations were assumed to be specific to each of the three measurement types (different procedure were used for their obtention).
They were assigned vague log-uniform prior distributions and were calibrated together with the other model parameters. A total of 11 parameters (8 kinetic parameters and 3 statistical ones) were calibrated. Posterior statistical distributions of the parameter values were obtained by Markov chain Monte Carlo (MCMC) sampling [5].
For the coupled PK-Nrf2 pathway model, the parameters directly controlling CsA kinetics were set to the joint posterior distribution mode found by the above calibration (see Table 1).
Another set of 27 structural parameters ( were computed the same way as the experimental ones, using the actual quantity predicted at the first time of the control experiment as a reference. The prior parameter distributions chosen were either vague or centered around the values used by Zhang [15] (see Table 2). The data likelihoods were assumed to be lognormal distributions around the model predictions.

Results
Results for the in vitro pharmacokinetic model have been previously reported in Wilmes et al.    As for the previous 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 (less than exponentially) while the Nrf2 protein quantity reaches systematically a plateau. Figure 6 shows The model 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 retrocontrol of ROS scavenging by ROS themselves, through the Nrf2 signaling pathway, induces a highly nonlinear behavior illustrated on Figure 6. ROS generation runs out of control at CsA

Supplementary Material
Differential equations of the Nrf2 model.

SUPPLEMENTARY MATERIAL
Differential equations of the Nrf2 model: