Estimation of the Mechanism of Adrenal Action of Endocrine-Disrupting Compounds Using a Computational Model of Adrenal Steroidogenesis in NCI-H295R Cells

Adrenal toxicity is one of the major concerns in drug development. To quantitatively understand the effect of endocrine-active compounds on adrenal steroidogenesis and to assess the human adrenal toxicity of novel pharmaceutical drugs, we developed a mathematical model of steroidogenesis in human adrenocortical carcinoma NCI-H295R cells. The model includes cellular proliferation, intracellular cholesterol translocation, diffusional transport of steroids, and metabolic pathways of adrenal steroidogenesis, which serially involve steroidogenic proteins and enzymes such as StAR, CYP11A1, CYP17A1, HSD3B2, CYP21A2, CYP11B1, CYP11B2, HSD17B3, and CYP19A1. It was reconstructed in an experimental dynamics of cholesterol and 14 steroids from an in vitro steroidogenesis assay using NCI-H295R cells. Results of dynamic sensitivity analysis suggested that HSD3B2 plays the most important role in the metabolic balance of adrenal steroidogenesis. Based on differential metabolic profiling of 12 steroid hormones and 11 adrenal toxic compounds, we could estimate which steroidogenic enzymes were affected in this mathematical model. In terms of adrenal steroidogenic inhibitors, the predicted action sites were approximately matched to reported target enzymes. Thus, our computer-aided system based on systems biological approach may be useful to understand the mechanism of action of endocrine-active compounds and to assess the human adrenal toxicity of novel pharmaceutical drugs.


Introduction
Because steroid hormones play an important role in a wide range of physiological processes, the potential to disturb endocrine effects is a major concern in the development of novel pharmaceutical drugs such as etomidate and aminoglutethimide [1]. The adrenal gland is the most common target for toxicity in the endocrine system in vivo, because steroid hormones are primarily synthesized through enzymatic reactions in the adrenal cortex [2][3][4][5]. Indeed, in these studies based on chemically induced endocrine lesions observed in vivo, the most frequent site of reported effects was the adrenal gland. Therefore, the prediction of human adrenal toxicity based on the mechanism of on-or off-target actions in the early stages of drug development is important.
The NCI-H295R human adrenocortical carcinoma cell line has been used to elucidate mechanisms of adrenal steroidogenic disrupting compounds [1,6]. The H295R cell line was established by Gazder and his collaborators in 1990 [7], which expresses all key steroidogenic enzymes and steroidogenesis-related proteins [7][8][9]. H295R cells have the physiological characteristics of zonally undifferentiated human fetal adrenal cells and the ability to produce steroid hormones found in the adult adrenal cortex [1,7,9]. In vitro 2 Journal of Toxicology bioassays using the H295R human cell line have been able to evaluate the effects of chemicals on steroid hormone production [10][11][12][13][14][15], steroidogenic enzyme activities [11,16,17], and the expression of steroidogenic genes [11,18]. In transcriptome studies, the mechanisms of action of many steroidogenic disrupting compounds have been qualitatively assessed in terms of adrenal toxicity. However, gene expression does not always reflect the production of steroid hormones [19]. Furthermore, measuring a few specific steroid hormones may not be a useful approach to study the mechanisms of steroidogenic disrupting effects in complex pathways such as adrenal steroidogenesis. To systematically understand how exogenous compounds affect adrenal steroidogenesis, simultaneous determination of all detectable steroid hormones and integrative analysis of these complex data would be important. As an exploratory approach to analyze complex data, ToxClust developed by Zhang and colleagues in 2009 is able to visualize concentration-dependent response relationships in the characteristics of chemically induced toxicological effects [20]. However, this exploratory approach is unable to provide a quantitative understanding of the mechanism of action of adrenal toxicants or reveal systematic information about the effect of each enzymatic reaction, interactions, and feedback in the adrenal steroidogenesis pathway.
Systems biology based on computational models of biological processes and the comprehensive measurement of biological molecules is the most powerful approach to quantitatively understand the influence of each factor in complex biological pathways. In recent studies by our collaborators, a computational model of adrenal steroidogenesis has been developed in NCI-H295R cells, including the steroidogenic disrupting effects of metyrapone to inhibit enzymatic reactions of CYP11B1 [21,22]. The model reproduces the dynamics of adrenal steroidogenesis in NCI-H295R cells and the influence of metyrapone. A current computational model of adrenal steroidogenesis was incorporated with a reaction of oxysterol synthesis as a bypass to consume cellular cholesterol [22]. In addition, all reactions in this model are described by a kinetic equation of the first-order reaction [22]. It is difficult to quantitatively evaluate the influence of each protein in the complicated system of adrenal steroidogenesis using the reported models, because it is simple and any biochemical and cellular biological information is not sufficient. For example, to clearly understand the cause of the change from the differentially dynamic patterns of steroid hormones, it is necessary to consider the substrate inhibition of steroidogenic enzyme because most of steroidogenic enzymes recognize multiple steroids as the enzymatic substrate. However, the substrate inhibition of steroidogenic enzyme cannot be described by the mathematical model based on kinetic equations of first-order reaction that does not consider Michaelis constant expressing the affinity of the substrate. To quantitatively estimate the mechanism of steroidogenic disrupting compounds from comprehensive experimental data of adrenal steroidogenesis in NCI-H295R cells, the reported model should be improved according to the following two points. First, the kinetic equation of enzymatic reactions should be exchanged from the first-order equation to a steady-state kinetic equation based on the mechanism of the enzymatic reaction. Because a mathematical model organized by first-order equations operates in a simple structuredependent manner, it does not show complex behavior based on molecular interactions, feedback, or regulation. Second, intracellular localization processes of cholesterol should be incorporated as a considerable mechanism. Because intracellular cholesterol molecules are stored as cholesterol esters or widely distributed as membrane components, only a few cholesterol molecules localized on the mitochondrial inner membrane are available for the adrenal steroidogenesis pathway [23,24]. Moreover, cholesterol-trafficking processes from the outer to inner mitochondrial membranes, which are regulated by steroidogenic acute regulatory (StAR) protein, are one of the rate-limiting steps in adrenal steroidogenesis [24]. By overcoming these limitations in the reported steroidogenesis model, systems analysis of adrenal steroidogenesis in H295R cells may be able to quantitatively estimate the mechanism of action of steroidogenic disrupting compounds.
In the present study, to quantitatively estimate the toxicological mechanism of endocrine-active compounds in adrenal steroidogenesis and to predict human adrenal toxicity of novel pharmaceutical drugs in the drug discovery phase, we developed a novel computational model of steroidogenesis in NCI-H295R cells. It includes cholesterol transport into intracellular regions from the extracellular space, the cholesterol translocation system in intracellular regions, including oxysterol synthesis, the metabolic pathway of adrenal steroidogenesis, and transport of steroid hormones. Global sensitivity analysis of this adrenal steroidogenesis model is able to evaluate the influence of each steroidogenic enzyme and related protein for each steroid hormone observed in an in vitro steroidogenesis assay of NCI-H295R cells. Furthermore, the mechanisms of action of steroidogenesis disrupting compounds for steroidogenic enzymes can be estimated by the optimization method to solve the reverse problem from the concentration changes of 12 steroid hormones measured by liquid chromatography/mass spectrometry in the steroidogenesis assay of NCI-H295R cells in vitro. Using this developed model of adrenal steroidogenesis and the analytical approach, the in vitro steroidogenesis assay of NCI-H295R cells can assess the human adrenal toxicity of a novel pharmaceutical drug based on quantitative understanding of its toxicological mechanism in adrenal steroidogenesis.

Adrenal Steroidogenesis in Human Adrenal
Corticocarcinoma NCI-H295R Cells. NCI-H295R cells were stimulated with adrenocorticotrophic hormone (ACTH), forskolin, and angiotensin II to initiate steroidogenesis. Changes in steroid concentrations over time were measured after stimulation in both cells and culture medium to construct a simulation model.

Mass Spectrometry.
A triple quadrupole mass spectrometer API4000 (Applied Biosystems/MDS Sciex, Concord, Canada) coupled with an electrospray ionization source was operated in the positive ion mode. The optimized ion source conditions were as follows: collision gas, 6 psi; curtain gas, 40 psi; ion source gas 1, 50 psi; ion source gas 2, 80 psi; ion source voltage, 5500 V; ion source temperature, 600 ∘ C. Nitrogen was used as the collision gas in the multiple reaction monitoring (MRM) mode. The conditions of declustering potential, collision energy, and collision cell exit potential were optimized by every steroid.

Estimation of the Cell Volume.
Cell volume was estimated from the number of cells in the well and the average diameter of the cells. Cells were detached from the well using 0.025% trypsin (MP Biomedicals, Inc., Irvine, CA) in a 0.02% EDTA solution (Dojindo Laboratories, Kumamoto, Japan) at the start of preculture, start of stimulation, and at 24, 48, and 72 h after stimulation. The numbers and diameters of the cells were measured by a cell counter Vi-cell XR 2.01 (Beckman Coulter, Krefeld, Germany) after trypan blue staining. Parameters of the cell volume and number of cells were estimated to fit experimental time-course data using exponential curves.

Validation Study Using Adrenal
Toxicants. NCI-H295R cells were cultured for 3 days in 6-well plates and then stimulated with the above-mentioned compounds. Upon the start of stimulation, various concentrations of test chemicals were added to the cultures. After a further 3 days of culture with the chemicals, the concentrations of 12 steroids (PREG, HPREG, DHEA, PROG, HPROG, DIONE, DCORTICO, DCORT, CORTICO, CORT, ALDO, and TESTO) in the culture medium were measured by LC/MS/MS. The test concentrations of the chemicals were determined by dose-finding cytotoxicity assays. The cytotoxicity assay was conducted in 96-well plates using ATP content in cells as an endpoint (CellTiter-Glo™ Luminescent Cell Viability Assay, Promega). Concentrations that caused more than 20% cytotoxicity were not used in the steroidogenesis assay. The test concentrations of adrenal steroidogenesis inhibitors and other compounds are shown in Table 1.

Statistical Analysis.
Comparisons were performed by the two-sample Welch's -test with Bonferroni multiple testing correction for each steroid hormone species. Statistically significant steroid hormones were considered at adjusted values of less than 0.01. Differential metabolic steroid profiles were classified by hierarchical cluster analysis. Pairwise distances between all compounds and all steroids were calculated by standardized Euclidean metric. This distance matrix was analyzed with Ward's method for hierarchical clustering. Statistical analysis was performed using MATLAB software (MathWorks, Inc., Natick, MA).

Mathematical Modeling of Adrenal Steroidogenesis in
NCI-H295R Cells. Steroid hormones secreted from human adrenal corticocarcinoma NCI-H295R cells are synthesized from cholesterol through the C 21 -steroid hormone biosynthesis pathway. A mathematical model of adrenal steroidogenesis in NCI-H295R cells was constructed with cholesterol transport and the intracellular localization pathway, the oxysterol synthesis pathway as a bypass of steroidogenesis, the C 21 -steroid hormone biosynthesis pathway as the main steroidogenesis pathway, passive transport of steroid hormones, and cell proliferation ( Figure 1). In this model, two compartments, the intracellular space and culture medium, were incorporated as the available region. Equations and parameters of the cell proliferation and diffusional transport of steroid hormones have been proposed by previous studies [21,22]. Cholesterol transport and the intracellular localization pathway including the oxysterol bypass were integrated using a part of the ACTH-stimulated cortisol secretion model described by Dempsher and colleagues [47]. The C 21 -steroid hormone biosynthesis pathway includes 14 steroid hormones, PREG, HPREG, DHEA, PROG, HPROG, DIONE, TESTO, DCORTICO, DCORT, CORTICO, CORT, ALDO, E1, and E2, and 17 enzymatic reactions catalyzed by nine steroidogenic enzymes, cholesterol side chain cleavage enzyme (CYP11A1), 17 -hydroxylase (CYP17H), C 17,20 -lyase (CYP17L), 3 -hydroxysteroid dehydrogenase (HSD3B2), 21hydroxylase (CYP21A2), 11 -hydroxylase (CYP11B1), 18hydroxylase (CYP11B2), 17 -hydroxysteroid dehydrogenase (HSD17B3), and aromatase (CYP19A1). In this mathematical model of adrenal steroidogenesis in NCI-H295R cells, the flux velocities of molecular transportation and enzymatic reaction rates of steroidogenic enzymes were defined based on the first-order reaction and rapid-equilibrium enzyme kinetics, respectively. All equations in the mathematical model of adrenal steroidogenesis of NCI-H295R cells were described in a supplementary document (see Supplementary Material available online at http://dx.doi.org/10.1155/2016/4041827). The rate constants and the maximum activities were estimated by fitting to experimental time-course data of the concentrations of cholesterol and all steroids. Initial values of cholesterol and the 14 steroid concentrations were used in each experimentally measured value, and every steroid concentration was assumed to rapidly reach the equilibrium state between the culture medium and intracellular space. All fixed values of static parameters and initial values of variable parameters in this model were described in Tables S1 and S2 in a supplementary document, respectively.

Modeling and Simulation Environment.
This computational model of adrenal steroidogenesis in NCI-H295R cells was developed on the simBio platform which is a general environment of biological dynamic simulation and computational model development [48]. ODEs were solved by the fourth-order Runge-Kutta method with a variable time step. The time step ( ) was adjusted to refer to the maximum absolute value of flux velocities or enzymatic reaction rates at each time point, and the range of the time step was from 1 × 10 −5 to 10 −2 . To confirm whether the range of the time step was suitable, the numerical error ratio was calculated by certain fixed time steps in the range of the time step, which was under 1 × 10 −8 in every time step. The duration time of computational simulation of adrenal steroidogenesis in NCI-H295R cells was set at 72 h.

Parameter Optimization.
To reconstruct experimental time-course patterns of the concentrations of cholesterol and the 14 steroids in the culture medium and intracellular space, we optimized every rate constant and maximum velocity of the steroidogenic enzymes. This parameter optimization problem was solved by the Levenberg-Marquardt method which is one of the nonlinear least squares methods [49][50][51]. The objective function of optimization was used as the following normalized least squares distance (NLSD): where ℎ is the compartment (culture medium or intracellular space), is the molecular species (cholesterol and the    in compartment ℎ at time point , sim ℎ, , is the simulated concentration of molecule in compartment ℎ at time point , and max ℎ, is the maximum concentration of molecule in compartment ℎ over all time points. Data points under the lower quantitation limit were excluded from the evaluation by the objective function.
Effects of every static model parameter for parameter optimization were calculated from differences of fitting the objective function using sensitivity analysis.

Quantitative Estimation of the Mechanism of Action of
Adrenal Toxicants. Metabolic steroid profiling and differential patterns of the adrenal steroid hormones by chemical perturbation were reconstructed to optimize the relative activities of the steroidogenic enzymes. The input data for the quantitative mechanistic analysis of adrenal toxic compounds was a fold change (ratio) of the measured 12 steroid concentrations induced by drug exposure for 72 h. The two-step optimization method of the real-coded genetic algorithm (RCGA) was adopted as a global optimization method in the quantitative mechanistic analysis of adrenal toxic compounds. The operations of the crossover and generation alteration model in RCGA were used for the real-coded ensemble crossover (REX) and just generation gap (JGG) [52][53][54][55]. As the initial parameters of RCGA, maximum generation, population size, selection size of parent individuals, population size of child individuals, and termination criteria were 1000, 100, 6, 25, and under 0.1 of NLSD, respectively. The search space for the relative activities of the steroidogenic enzymes was from 1/100 to 100. To evaluate the fitness of each individual, the sum of squared residuals for fold changes of measured 12 steroid concentrations was used as the objective function. Nonlinear least squares optimization by the Levenberg-Marquardt method was used as a local search [49][50][51]. As the estimated mechanisms of actions of the adrenal toxic compounds, the relative activities of eight steroidogenic enzymes (CYP11A1, CYP17H, CYP17L, HSD3B2, CYP21A2, CYP11B1, CYP11B2, and HSD17B3) were optimized by the above-mentioned 2-step optimization method. Every optimization calculation was duplicated to check the numerical stability of the optimal parameters.

Global Dynamic Sensitivity
Analysis. The property of every kinetic parameter in this computational model of steroidogenesis in NCI-H295R cells was evaluated by dynamic sensitivity analysis. The sensitivity ( , ) of kinetic parameter for variable parameter was defined by the following equation: where variable parameter was the concentration of a steroid hormone in the cytosolic space of NCI-H295R cells. The perturbation for kinetic parameters was +10% (Δ / = 0.1).

Differentially Steroid Profiling of Adrenal Toxicants.
After NCI-H295R cells were exposed to each test compound during three days, the concentrations of 12 steroid hormones in the culture medium were simultaneously measured by LC/MS/MS. All effects of the compounds on adrenal steroidogenesis were evaluated at the concentration without any overt cytotoxicity. The differential metabolic steroid profiles of 11 adrenal toxic compounds were classified and visualized by using hierarchical clustering analysis ( Figure 3).
Four adrenal toxicants without steroidogenic inhibition, AN, SM, TG, and FN, did not change the medium concentrations of all steroid hormones by more than 2-fold. Abovementioned 4 compounds at every condition and 7 adrenal steroidogenic inhibitors at the low exposure concentration were gathered into a big cluster as nonchange group. The 7 steroidogenic inhibitors at the maximum exposure concentration showed the characteristic steroid profiles each, but 100 M DZ and 10 M SP were classified as a cluster.

Optimization of the Mathematical Model of Adrenal
Steroidogenesis in NCI-H295R Cells. The mathematical model of adrenal steroidogenesis in NCI-H295R cells was optimized for several kinetic parameters of cholesterol transport, intracellular localization, the oxysterol pathway, and maximum velocity of steroidogenic enzymes to fit the experimental time-course data. All optimized kinetic parameters are shown in Table S1 in a supplementary document. The optimized mathematical model was reconstructed with the experimental dynamic patterns of cholesterol and the 14 steroid hormones in the intracellular space and culture medium. The fitness was 0.621761 of NLSD values as the fitting objective function. The simulation results and experimental data are shown in Figure 4.
Optimized kinetic parameters were calculated sensitivities for the NLSD value as the fitting score and are shown in Table S1 in a supplementary document. The highly sensitive parameters for fitting the NLSD score were the extracted nine kinetic parameters, Cholesterol Transport , acc , loc , loc , CYP11A1 max , CYP11A1 , CYP17H maxA , HSD3B2 maxA , and CYP21A2 maxA , which had higher than 3.0 fitting sensitivity.

AGT.
The mechanism of action of AGT in adrenal steroidogenesis was estimated by inhibition of CYP11A1, HSD3B2, CYP21A2, and CYP11B1 at 100 M (estimated inhibitions were 77.0%, 78.0%, 81.1%, and 59.8%, resp.) ( Figure 5(e)). AGT has been reported to inhibit CYP11A1, CYP21A2, CYP11B1, and CYP11B2 [6,[27][28][29][30]. Our results were mostly consistent with the previous reports. In particular, CYP11A1 appeared to be inhibited strongly by AGT. In our study, HSD3B2 inhibition of AGT was shown by mechanistic analysis based on systems biology approaches as a novel mechanism of action of AGT. Inhibition of AGT by CYP11B2 was not estimated in our study. However, the concentration of ALDO in the culture medium decreased to 3.8% of the normal stimulated condition. Inhibition of AGT by CYP11B2 has been shown using sheep adrenal homogenates as well as a human adrenal homogenate from a patient with Cushing's syndrome [30]. The activity of 18-hydroxylase induced by CYP11B2 was determined as the conversion of corticosterone to 18-hydroxycorticosterone in the previous study. The cause of the discrepancy regarding the effect of AGT on CYP11B2 was suggested to be substrate inhibition, because the intracellular concentration of CORTICO was increased by over 10 times of that in the culture medium to reach 50 M. Another possibility was poor quantitative reliability of the experimental data, because the ALDO concentration was under the lower limit of quantification at 100 M AGT. Hecker and colleagues reported that 3 M AGT decreases PREG and PROG concentrations and increases the TESTO concentration [10]. However, AGT did not increase the TESTO concentration in our study. One possibility is that the concentration of TESTO was already enhanced by about 3.3-fold through stimulation with ACTH, forskolin, and angiotensin II.

DDD.
The mechanism of action of DDD in adrenal steroidogenesis was estimated by dose-dependent inhibition of CYP11A1, HSD3B2, CYP21A2, and CYP11B1 (estimated inhibitions at 25 M were 87.0%, 86.9%, 76.9%, and 84.9%, resp.) ( Figure 5(f)). DDD has been reported to inhibit CYP11A1, HSD3B2, CYP21A2, CYP11B1, and CYP11B2 [29,31,32]. Inhibition of DDD by CYP11B2 was not estimated in our study. However, the concentration of ALDO in the culture medium decreased to 3% of that in the normal stimulated condition. Inhibition of DDD by CYP11B2 has been shown using mitochondrial and microsomal fractions prepared by standard centrifugation procedures from a bovine adrenal cortex homogenate [32]. The cause of the discrepancy regarding the inhibition of DDD by CYP11B2 could not be explained by same effect in the case of AGT.

SP.
The mechanism of action of SP in adrenal steroidogenesis was estimated by inhibition of HSD3B2, CYP21A2, and HSD17B3 (estimated inhibitions at 10 M were 70.2%, 59.5%, and 59.3%, resp.) ( Figure 5(g)). SP has been reported to inhibit CYP17H, CYP17L, CYP11B1, and CYP11B2 [6,[33][34][35]. The inhibitory effect of SP on the HSD3B2 enzyme was a novel mechanism of action. The main action of SP is as a mineralocorticoid receptor (MR) antagonist. SP has also been reported to exert some off-target effects by binding to androgen, glucocorticoid, and progesterone receptors [56][57][58]. SP has been shown to inhibit the production of ALDO and CORT from PREG induced by angiotensin II in H295R cells, but the specific MR antagonist eplerenone did not show the inhibitory effects [59]. Therefore, HSD3B2 inhibition by SP is not mediated via MR, and the action might be direct inhibition of HSD3B2 enzymes or a part of known off-target effects mediated through other nuclear hormone receptors. Regarding CYP17H and CYP17L, our results were consistent with previous reports [33,34]. 7 -thiospironolactone, which is synthesized by deacetylation of SP, inhibits CYP17H and CYP17L [34]. The fact that there were no inhibitions of CYP17H or CYP17L in our study suggests that SP might not be deacetylated to 7 -thiospironolactone in NCI-H295R cells.
Regarding CYP11B1 and CYP11B2, our results were unclear compared with a previous study. It has been shown that 30 M SP inhibits CYP11B1 and CYP11B2 in human and bovine adrenal mitochondria [35]. The cause of CYP11B1 and CYP11B2 inhibition by SP could not be determined in our study, which might be due to the lower maximum exposure concentration of SP than that in the previous report. We could not examine SP concentrations over 10 M because these concentrations were cytotoxic in NCI-H295R cells.

MP.
The mechanism of action of MP in adrenal steroidogenesis was estimated by dose-dependent inhibition of CYP11B1 (estimated inhibitions at 1, 10, and 100 M were 57.1%, 82.7%, and 98.2%, resp.) ( Figure 5(h)). MP has been reported to inhibit CYP11B1 as its major effect and CYP11A1 and CYP11B2 as a weak effect [6,29,[36][37][38][39]. The results were able to show that MP is a selective inhibitor of CYP11B1 in the previous report. However, the estimated effect of MP at a high concentration, 100 M as the maximum exposure concentration, was unclear. According to the previous report, selectivity of MP for CYP11B1/CYP11B2 is about five times [39]. In addition, 20 M MP has a slight inhibition effect on CYP11A1 in H295R cells [29].

MC.
The mechanism of action of MC in adrenal steroidogenesis was estimated by inhibition of CYP17H, CYP17L, CYP11B1, and HSD17B3 (estimated inhibitions at 10 M were 69.1%, 53.0%, 76.4%, and 57.1%, resp.) ( Figure 5(j)). MC has been reported to inhibit not only CYP17H and CYP17L but also CYP11A1, CYP21A2, and CYP11B1 [41,44,45]. The results in the previous reports were able to estimate that MC is a CYP17 inhibitor. However, CYP11A1 inhibition by MC, probably instead of a reduction in StAR expression, was not clearly detected in our study using NCI-H295R cells, because there were no decreases in the concentrations of PREG and PROG in the culture medium. Indirect inhibition of CYP11A1 via the peripheral-type benzodiazepine receptor has been reported in mouse adrenocortical Y-1 cells treated with MC in the absence of stimuli by measuring PREG production [44]. On the other hand, reductions of StAR protein expression and/or transport activity without affecting total steroid synthesis or CYP11A1 and HSD3B2 enzyme expression or activities have been reported in (BU) 2 cAMPstimulated MA-10 Leydig tumor cells treated with MC by measuring PROG production [45]. Therefore, the effect of MC on the initial reaction in adrenal steroidogenesis from cholesterol should be different according to the cell type and stimulation condition. Inhibition of CYP21A2 and CYP11B1 by MC has been reported as decreases in the consumption of PROG and DCORTICO, respectively [41]. Inhibition of CYP11B1 was estimated by the action of MC in this study, but that of CYP21A2 was not detected. In the previous experimental report, inhibitory sites by MC might have been reflected by inhibition of CYP17H activity, because CYP21A2 activity was measured as a decrease in labeled PROG.

DZ.
The mechanism of action of DZ in adrenal steroidogenesis was estimated by inhibition of CYP11A1, HSD3B2, CYP21A2, CYP11B1, and HSD17B3 (estimated inhibitions at 100 M were 58.6%, 94.1%, 96.5%, 87.2%, and 98.1%, resp.) ( Figure 5(k)). DZ has been reported to inhibit HSD3B2 and CYP21A2 [46]. The results of HSD3B2 and CYP21A2 were consistent with the previous report. However, inhibitions have not been reported for CYP11A1, CYP11B1, and HSD17B3. These estimated effects of DZ on CYP11B1 and HSD17B3 were unclear, because the concentrations of ALDO and TESTO were less than the lower limit of quantification at 100 M. In addition, these enzymes act downstream of the strong action points of DZ, such as HSD3B2 and CYP21A2.

Dynamic Sensitivity Analysis of Adrenal Steroidogenesis.
To comprehensively understand the dynamics of adrenal steroidogenesis, dynamic sensitivities were calculated for steroid concentrations secreted by NCI-H295R cells using our constructed mathematical model of steroidogenesis. The results of dynamic sensitivity analysis at 72 h of duration and 6 h of interval time are presented as a heat-map in Figure 6. The top 10 parameters of the total area under the curve of dynamic sensitivity for cholesterol and the 14 steroids in culture medium were HSD3B2 maxA , CYP11A1 max , CYP21A2 maxA , loc , loc , CYP11A1 , CYP17H maxA , Cholesterol Transport , acc , and CYP21A2 maxB in order from the top. Cholesterol uptake ( Cholesterol Transport ), StAR protein ( loc ), and CYP11A1 ( CYP11A1 max ), which are determining factors of the capacity for steroidogenesis, promoted the production of mineralocorticoids (DCORTICO, CORTICO, and ALDO) and restrained the synthesis of glucocorticoids (DCORT and CORT) and sex steroids (DIONE, TESTO, and E1) because of the accumulation of intermediate molecules in steroidogenesis (PREG, HPREG, DHEA, PROG, and HPROG) only by self-activation. The dynamic patterns of the intermediate molecules in steroidogenesis were mainly dependent on the activity of CYP17H and HSD3B2 with PREG as the substrate of these enzymes, in which the dynamic sensitivities of CYP17H maxA for HPREG, and HPROG and HSD3B2 maxA for PROG, HPROG, and DCORTICO reversed the direction of sensitivity at 49-66 h after stimulation. The dynamic sensitivities of the maximum activities of HSD3B2 for PREG ( HSD3B2 maxA ) and CYP21A2 for PROG ( CYP21A2 maxA ) were related to all steroids at 72 h. Almost all model parameters had positive sensitivity for downstream steroids in the adrenal steroidogenic pathway and negative sensitivity for direct-binding steroids as substrates of the steroidogenic enzyme. The sensitivity of max in all steroidogenic enzymes was relatively higher than for the same steroid substrate.

Simulation of the Metabolic Balance of Adrenal
Steroidogenesis Pathway. To clearly show the property of the metabolic shift between mineralocorticoid and glucocorticoid biosynthesis, we performed two-dimensional parameter scanning of the enzymatic activities of CYP17H and HSD3B2 (Figure 7). NCI-H295R cells lost the ability to produce all steroid hormones when enzymatic activities of CYP17H and HSD3B2 were changed by over 60% and 30%, respectively. Activation of CYP17H and/or HSD3B2 induced the metabolic shift that enhanced the glucocorticoid biosynthesis and deviated from the mineralocorticoid biosynthesis. On the other hand, inhibition of CYP17H and/or HSD3B2 induced the metabolic shift that enhanced the mineralocorticoid biosynthesis and deviated from the glucocorticoid biosynthesis. Moreover, the enzymatic activity of HSD3B2 regulated the metabolic balance of sex steroids and the precursors on adrenal steroidogenesis of NCI-H295R cells. E1, TESTO, and DIONE were produced by NCI-H295R cells when activating the enzymatic activity of HSD3B2. Conversely, E2 and DHEA were produced by NCI-H295R cells when suppressing the enzymatic activity of HSD3B2. The biosynthesis of downstream steroids in adrenal steroidogenesis pathway, such as mineralocorticoids and glucocorticoids, was almost completely terminated when the enzymatic activity of HSD3B2 was decreased by over 80%.

Importance of 3 -HSD Activity in Adrenal Steroidogenesis.
Our systematic analysis using the mathematical model of adrenal steroidogenesis in NCI-H295R cells revealed that the enzymatic activity of 3 -HSD controls the dynamics of adrenal steroidogenesis. The activity of the StAR protein controls the net capacity of steroidogenesis in steroidogenic cells, which is the transport of cholesterol from the outer to inner mitochondrial membranes. Both the expression levels of StAR protein and mRNA are rapidly elevated in response to stimulation by tropic hormones such as ACTH [62,63]. Another important factor in adrenal steroidogenesisis is the cholesterol side chain cleavage enzyme CYP11A1, the first rate-limiting and hormonally regulated step in the synthesis of all steroid hormones, which is conversion of cholesterol to pregnenolone in mitochondria [64]. According to our results of global sensitivity analysis (Supplementary Table 1 and Figure 6(d)), in addition to CYP11A1 and StAR proteins, 3 -HSD was one of the key regulators in adrenal steroidogenesis of NCI-H295R cells. And also, this result suggests that a significant regulatory mechanism in steroidogenesis 14 Journal of Toxicology pathway is very reasonable. StAR, CYP11A1, and 3 -HSD (isoforms 1 or 2 in humans) proteins generally respond to the same hormones that stimulate steroid production through common pathways such as cAMP signaling in adrenal glands and testes [65,66]. Moreover, our data also support recent experimental evidence from clinical and in vivo studies, suggesting that the enzymatic activity of 3 -HSD plays an important role in the regulation of mineralocorticoid synthesis in adrenal steroidogenesis and contributes to hypertension caused by abnormal overproduction of aldosterone [67][68][69][70]. Circadian clock-deficient Cry-null mice show saltsensitive hypertension due to abnormally high synthesis of aldosterone, which is caused by constitutively high expression of HSD3B6 mRNA and protein in the adrenal cortex [67,68]. Recent clinical observations have revealed predominant expression of HSD3B2 mRNA and protein in tumor cells of aldosterone-producing adenoma (APA), and HSD3B1 mRNA significantly correlated with CYP11B2 mRNA levels and plasma aldosterone concentrations in APA patients [69,70]. However, the relationship is unclear and disputed in a small-scale clinical study indicating that genetic variation in HSD3B1 affects blood pressure and hypertension in APA patients [71]. The results of our simulation study suggest that 3 -HSD protein (human genes are HSD3B1 and HSD3B2) is one of the determination factors for the dynamic property of adrenal steroidogenesis. Our results support the clinical evidence of Doi and colleagues [69], and we believe that the HSD3B1 enzyme has a promising potential as novel drug target for endocrine hypertension.

Metabolic Shift of Adrenal
Steroidogenesis and the Contributions of HSD3B2 and CYP17H. The metabolic properties of adrenal steroidogenesis in NCI-H295R cells were revealed by dynamic sensitivity analysis using the mathematical model ( Figure 6). Mineralocorticoids, such as DCORTICO, COR-TICO, and ALDO, and intermediate steroids upstream of the adrenal steroidogenesis pathway, such as PREG, HPREG, DHEA, PROG, and HPROG, were accelerated by reactions of cholesterol import ( Cholesterol Transport ), StAR protein ( MTR and loc ), and CYP11A1 ( CYP11A1 max ). On the other hand, glucocorticoids, such as DCORT and CORT, and sex hormones, such as DIONE, TESTO, and E1, were suppressed by these model parameters. Therefore, enhancement of the net adrenal steroidogenesis capacity, which supplies PREG precursor to the pathway, causes a production shift from glucocorticoids to mineralocorticoids by substrate inhibitions of CYP17H, HSD3B2, and CYP21A2 caused by accumulation of initial intermediate steroids such as PREG and PROG. Sensitivities of CYP17H ( CYP17H maxA ) and HSD3B2 ( HSD3B2 maxA ) for the products were dynamically changed and these parameters determined the metabolic balance of downstream steroids in the adrenal steroidogenesis pathway. According to these results of dynamic sensitivity analysis of StAR, CYP11A1, CYP17H, and HSD3B2, we suggest that the enhancement of CYP17H and HSD3B2 activity during ACTH stimulation was important to shift the steroidogenic output away from ALDO biosynthesis towards CORT biosynthesis, as well as adrenal androgen production. This suggestion partially Journal of Toxicology supports a comparative animal study in which molecular and cellular variations in CYP17H activity dramatically affect acute cortisol production, resulting in distinct physiological and behavioral responses [72].
Results of two-dimensional parameter scanning of the enzymatic activities of CYP17H and HSD3B2 quantitatively showed the detail of the metabolic relationship between mineralocorticoid and glucocorticoid biosynthesis (Figure 7). Particularly, the results showed that the balance of these enzymatic activities was very important for the typical function of NCI-H295R cells, namely, the ability to produce all steroid hormones. NCI-H295R cells lost this function when enzymatic activities of CYP17H and HSD3B2 were changed by over 60% and 30%, respectively. In addition, they became mineralocorticoid (ALDO) secreting cells when the enzymatic activity of CYP17H or HSD3B2 was inhibited by over 50% or glucocorticoid (DCORT and CORT) secreting cells when these enzymes were activated by over 50%. In particular, this analysis also showed that HSD3B2 was a key player in the adrenal steroidogenesis of NCI-H295R cells, because HSD3B2 inhibition by over 80% almost completely inhibited the biosynthesis of downstream steroids. The ratio of CYP17A1 to HSD3B2 mRNA expression levels has been related to several endocrine diseases with a low level in APAs [73] and high level in cortisol-producing adenomas [74]. Furthermore, the expression levels or enzymatic activities of CYP17A1 and HSD3B1 have been related to androgen production in polycystic ovary syndrome [75,76]. These clinical studies support our simulation results indicating that the balance of enzymatic activity of CYP17H and HSD3B2 determines the shift in steroidogenic output to mineralocorticoids, glucocorticoids, or androgens.

Methodologies of Quantitative Mechanistic Analysis for
Drug Discovery. According to our results obtained using the mathematical model of steroidogenesis in NCI-H295R cells, such as sensitivity analysis, comprehensive analysis based on systems biology is available to quantitatively estimate the mechanism of action of steroidogenic disrupting compounds from differential profiling of adrenal steroid hormones, because dynamic patterns of steroid hormones in adrenal steroidogenesis pathway are highly complex. Our proposed method of quantitative mechanistic analysis of steroidogenic inhibitors was able to predict known action sites in the adrenal steroidogenesis pathway at only one time point (72 h after drug exposure). Moreover, according to the results of sensitivity analysis (Figure 6), max of all steroidogenic enzymes was more sensitive than the , because the intracellular concentrations of steroid hormones were almost maintained at sufficiently high levels compared with values of steroidogenic enzymes. These results suggested that estimation of the mechanism of action of drugs is more effective and detectable when using the influences of max as the searching parameters such as our proposed method. Our data showed that the proposed method based on a systems biology model is a very powerful tool for exploratory screening of steroidogenic disrupting compounds.
RCGA as a solver of parameter estimation problems in systems biology has been applied to biological network identification of gene regulatory networks and metabolic pathways and optimization of biological processes using experimentally observed time-course data [77][78][79][80][81][82]. In this study, RCGA was useful to estimate the mechanism of action of novel pharmaceutical drug candidates for adrenal steroidogenesis as a new application of RCGA in systems biology. We had two issues when applying RCGA to the quantitative mechanistic analysis of drug actions. These issues were the vast calculation cost and multimodality of quasi-optimum solutions in solving the optimization problem, because the mathematical model in systems biology consists of many equations and parameters. A proposed optimization strategy using RCGA based on REX/JGG was a highly stable and efficient calculation method for a better quasi-optimum solution than the unimodal normal distribution crossover (UNDX)/minimum generation gap (MGG) method that is well applied in the engineering field. In addition, we expanded the RCGA optimization program based on REX/JGG to a hybrid method of GA and then applied a local search as recommended by Harada and Kobayashi [28,83]. A final optimal solution was obtained with a good convergence property. Because these problems are general in systems biology studies, we suggest that the proposed hybrid method based on REX/JGG is a very useful tool for quantitative mechanistic analysis of novel pharmaceutical drugs, not limited to steroidogenic disrupting compounds.

Conclusions
The novel mathematical model of adrenal steroidogenesis was constructed in this study, including cholesterol transport and distribution, the C 21 -steroid hormone pathway, steroid transport, and cell proliferation, which could reproduce adrenal steroidogenesis in NCI-H295R cells. According to the results of dynamic sensitivity analysis using the new model, HSD3B2 plays the most important role in the metabolic balance of adrenal steroidogenesis in NCI-H295R cells. Moreover, to quantitatively estimate mechanisms of action of adrenal toxic compounds, we analyzed differential metabolic profiles of 12 steroid hormones at 3 days after exposure to 11 adrenal toxic compounds, by using the new mathematical model and a hybrid optimization method of the RCGA and a local search (nonlinear least squares). We could estimate which steroidogenic enzymes were affected by these compounds using the hybrid optimization method. Vasculotoxic agents were estimated to have no effect according to the results obtained by our method. In terms of adrenal steroidogenic inhibitors, the predicted action sites were approximately matched to the target enzymes as reported in the literature. Thus, our computer-aided method based on a systems biology approach may be useful to analyze the mechanism of action of endocrine-disrupting compounds and to assess the human adrenal toxicity of novel pharmaceutical drugs based on steroid hormone profiling.