A Semi-Physiologically Based Pharmacokinetic Pharmacodynamic Model for Glycyrrhizin-Induced Pseudoaldosteronism and Prediction of the Dose Limit Causing Hypokalemia in a Virtual Elderly Population

Glycyrrhizin (GL) is a widely used food additive which can cause severe pseudoaldosteronism at high doses or after a long period of consumption. The aim of the present study was to develop a physiologically based pharmacokinetic (PBPK) pharmacodynamic (PD) model for GL-induced pseudoaldosteronism to improve the safe use of GL. Since the major metabolite of GL, glycyrrhetic acid (GA), is largely responsible for pseudoaldosteronism via inhibition of the kidney enzyme 11β-hydroxysteroiddehydrogenase 2 (11β-HSD 2), a semi-PBPK model was first developed in rat to predict the systemic pharmacokinetics of and the kidney exposure to GA. A human PBPK model was then developed using parameters either from the rat model or from in vitro studies in combination with essential scaling factors. Kidney exposure to GA was further linked to an Imax model in the 11β-HSD 2 module of the PD model to predict the urinary excretion of cortisol and cortisone. Subsequently, activation of the mineralocorticoid receptor in the renin-angiotensin-aldosterone-electrolyte system was associated with an increased cortisol level. Experimental data for various scenarios were used to optimize and validate the model which was finally able to predict the plasma levels of angiotensin II, aldosterone, potassium and sodium. The Monte Carlo method was applied to predict the probability distribution of the individual dose limits of GL causing pseudoaldosteronism in the elderly according to the distribution of sensitivity factors using serum potassium as an indicator. The critical value of the dose limit was found to be 101 mg with a probability of 3.07%.


Introduction
Glycyrrhizin (GL), one of the major components in licorice root, is widely used as a sweetener in food products and chewing tobacco. In addition, GL is of clinical importance as an antiinflammatory and hepatoprotective agent in the treatment of dermatitis, gastritis, hepatitis and other conditions [1]. However, several studies have reported that excessive or long-term consumption of GL can induce adverse effects such as hypertension, hypokalemia and sodium retention, all of which are symptoms of pseudoaldosteronism. In serious cases, the outcome can be death [2].
The mechanism of this pseudoaldosteronism seems to be clear [3]. Glycyrrhetic acid (GA), the main metabolite of GL, has been verified to be the major substance leading to pseudoaldosteronism. GA is a potent inhibitor of 11b-hydroxysteroiddehydrogenase 2 (11b-HSD 2), the enzyme that converts cortisol to cortisone in the kidney. This inhibition leads to an increase in cortisol levels as well as in the cortisol:cortisone ratio in plasma and urine and, as cortisol and aldosterone have the same affinity for the mineralocorticoid receptor (MR) [4], to increased stimulation of the MR leading to increased electrolyte levels and the symptoms of pseudoaldosteronism.
The metabolism of GL depends on its route of administration [5]. GL administered by the intravenous (i.v.) route is partially hydrolyzed in the liver by lysosomal b-D-glucuronidase to 3-monoglucuronyl-glycyrrhretinic acid (3MGA). GL and 3MGA can then be excreted into the bile via canalicular transporters and finally reach the intestinal lumen. There GL and 3MGA are metabolized by the intestinal microbiota to GA, which then undergoes enterohepatic recycling. When GL is administered by the oral route (p.o.), the intestinal microbiota bring about its hydrolysis to GA which again undergoes enterohepatic recycling.
In assessing the risk of a food additive, the traditional method of measuring the acceptable daily intake, following the rule of 'dividing by ten', is arbitrary and provides an inaccurate estimation. It is a feasible and recommended method when the mechanism of action for the compound in question is not quite clear. However, a large number of in vivo and in vitro studies have recently reported the mechanisms of GL pharmacokinetics and GL-induced pseudoaldosteronism. We therefore attempted to use the method of mechanism-based modeling to make the risk control process more accurate and reasonable. Moreover, according to the opinions of the Scientific Committee on Food (SCF) [6] and the Dutch Nutrition Information Bureau (DNIB) [7], the upper safe limits for regular ingestion of GL are 100 mg/day and 200 mg/day respectively. However, because human toxicity studies remain incomplete, these upper limits may actually be too high to ensure sufficient protection to some subgroups. Consequently, modeling and simulation of scenarios difficult to undertake in clinical studies are needed to provide more accurate predictions of the dose limit in various populations. In addition, since the kidney is the target tissue of GA in GL-induced pseudoaldosteronism, determination of kidney exposure to GA in such populations is important to predict the likelihood of adverse effects of GL.
Physiologically-based pharmacokinetic (PBPK) modeling is now an established technique for estimating the absorption, distribution, metabolism and excretion of a drug and, in combination with pharmacodynamic (PD) modeling, predicting its pharmacological and/or adverse effects. Ploeger and colleagues developed a PBPK/PD model of GL in rat and human [8,9,10] and successfully predicted GL and GA plasma concentrations and the urinary cortisol:cortisone ratio in human. However, in the Ploeger model, the urinary cortisol:cortisone ratio did not indicate whether an adverse event had actually occurred. We therefore recognized the need to include in the model the clinical indicators of pseudoaldosteronism (e.g. serum potassium and sodium levels) and the kidney exposure to GA in order to better predict the adverse effects of GL. Moreover, once the MR is overstimulated, the renin-angiotensin-aldosterone system (RAAS) is inhibited as a compensatory effect meaning RAAS also plays an important role in changing the levels of body fluids and electrolytes when GL is ingested.
In the present study, we first developed a semi-PBPK model for GL and GA in rat and then, based on parameters in the rat model and species scaling, developed a human PBPK model in which kidney exposure to GA was further combined with a biological system-PD model to produce a PBPK/PD model capable of predicting the indicators of GL-induced pseudoaldosteronism in human.

Experimental studies
During the modeling process, optimization and validation were conducted step by step using in vivo experimental datasets largely based on the literature. Moreover, in order to ask the validation more robust, we found it necessary to have data relating to the PK of GA and its kidney distribution after p.o. administration of GL. Two separate studies as follows were first conducted to provide a validation dataset for the rat PBPK model.

Materials
GL and GA (purity .98%) were obtained from Chia-tai Tianqing Pharmaceutical Co, Nanjing, China. Diclofenac for use as internal standard (IS) was purchased from Nanjing Chemical Reagent Co Ltd., Nanjing, China. All other chemicals were of analytical grade and used as received. Ultrapure water was used in all studies.
Male Sprague-Dawley rats (weight 220-250 g) were obtained from B&K Universal Group Ltd., Shanghai, China. The rats were housed under a 12 h light/ dark cycle with free access to commercial rat chow and water. All animals were fasted overnight with free access to water before an experiment. All animal experiments were approved by the Animal Ethics Committee of China Pharmaceutical University and performed under a license granted by the Jiangsu Science and Technology Office (China). All efforts were made to minimize suffering.
Pharmacokinetics of GA after p.o. administration of GL GL was dissolved in ultrapure water and administered p.o. to 5 rats at 100 mg/kg. Blood samples (0.2 ml) were collected at the following time points: 0.5, 1, 2, 4, 6, 8, 12, 24 h after the dose. Samples were stored at 220˚C pending analysis.
Kidney distribution of GA after p.o. administration of GL A group of 12 rats received a p.o. dose of GL at 200 mg/kg. After 9, 16 and 22 h, 4 animals were sacrificed at each time point by femoral artery exsanguination and blood and kidneys harvested. Each tissue was homogenized with an equal volume of cold saline after which samples were stored at 220˚C pending analysis.
Analysis of GA in rat plasma and kidney GA concentrations in the above samples were analyzed according to a previously developed method [11]. 90 ml plasma or tissue homogenate was taken into an Eppendorf tube, and added to 10 ml diclofenac (internal standard, 1 mg/ml) and 200 ml acetonitrile. The tube was vortex mixed and then centrifuged at 130006g for 10 min. The supernatant was transferred to another clean Eppendorf tube and centrifuged at 130006g for 10 min. An aliquot (5 ml) of supernatant was injected directly onto the HPLC column. And the analysis was carried out by LC-MS using an LC-MS 2010 system (Shimadzu, Kyoto Prefecture, Japan) comprising 2 LC-20AD pumps, a CBM-20A system controller, a CTO-20A column oven set at 40˚C, a SIL-20AC autosampler set to inject 5 mL and a single quadrupole mass spectrometer equipped with an electrospray ionization (ESI) interface operated in the negative ionization mode. Separation was carried out on an Agilent Microsorb 100-5 C18 column (15062.0 mm) using gradient elution with 0.1% acetic acid-10 mM ammonium acetate in water as mobile phase A and acetonitrile as mobile phase B at a fl ow rate of 0.2 mL/min. The gradient (min, %B) was as follows: 0-0.5, 20; 0.5-2.5, 40; 2.5-6.5, 90; 6.5-10.0, 20. Elution in the period 2.5-7.0 min was delivered into the mass analyzer. Detection was by selected ion monitoring (SIM) of ions with 469.3 for GA and 294.1 for diclofenac (internal standard IS). Optimum MS parameters were: Capillary voltage 2 1.5 kV; nebulizer nitrogen flow rate 1.5 L/min; drying gas temperature 280˚C; block temperature 220˚C.

Model building workflow
Following the proposed workflow shown in Figure 1, a PBPK/PD model for GLinduced pseudoaldosteronism was developed and used to predict the dose limit of GL producing pseudoaldosteronism in a virtual population. In brief, the rat PBPK model was first developed. Physicochemical parameters, absorption, distribution, metabolism and elimination (ADME) data, and physiological parameters were collected from the literature and used as initial input parameters of the model. After that, in vivo data sets were used: some to optimize, some to validate the model. If the model adequately reflected the observed data, the modeling process was continued; if not, the process would revert to the first step (see flowchart in figure 1). Based on the structure and parameters of rat PBPK model and also on scaling-induced extrapolation, we then developed the human PBPK model following similar steps as previously mentioned. After validation of the PBPK model, kidney exposure to GA in human was output and linked with the human PD model (consisting of three modules: 11b-HSD2 module, RAAS module and Electrolyte module). Finally, clinical indicators of pseudoaldosteronism could be predicted. The whole PBPK/PD model is developed step by step and validated during each step by observed data to prevent systematic errors.

Description of the GL model
The model consists of four tissue compartments and one blood compartment as shown in Figure 2A. The abbreviation of the parameters used in the PBPK models  Table 1. The kidney is the main target organ of GL-induced pseudoaldosteronism and has a high blood perfusion rate. The gut is isolated to describe the progress of absorption. Both of these tissues are assumed to be perfusion limited. As mentioned above, hepatic uptake and biliary excretion play important roles in GL elimination. GL has been shown to be actively transported into the liver by organic anion transport proteins (OATPs) [12] and subsequently excreted into the bile by multidrug resistance-associated protein 2 (MRP2) [13]. In order to model these processes, the liver was divided into a venous compartment and a tissue compartment and active transport of GL was described by a Michaelis-Menten-type equation. Besides biliary excretion, GL also undergoes passive diffusion from the liver into the plasma or is hydrolyzed to 3MGA by b-glucuronidase, both of which were assumed to be first order processes [14]. Then distribution in the remaining tissue was described by the parameters K 12_r and K 21_r . Since urinary excretion of GL and its metabolites accounts for only a small proportion (,2.5%) of the total drug [15,16], renal clearance was ignored for GL and its metabolites. Sum of small intestine and large intestine. c from literature [17]. d from literature [8,9].
doi:10.1371/journal.pone.0114049.t002 Table 3. Biochemical parameters for the physiologically based model for GL, GA and GAM in rat. Calibration of the GL model Anatomical and physiological parameters taken from the literature [17] are listed in Table 2. Partition coefficients for each compartment [18] and the unbound fraction f u of GL in plasma [19] are listed in Table 3. The unbound fraction of GL in tissue fu t was estimated by equation (1) where fu p is the unbound fraction in plasma, R is the average value of the tissue interstitial fluid-to-plasma ratio of albumin and lipoproteins (R50.5 for lean tissues and R50.15 for adipose tissue) [20,21]:

Parameters
The initial values of V max and K m for hepatic uptake were taken as the values determined in an isolated rat liver perfusion study [22]. As V max and K m can change simultaneously to keep the ratio constant, we here fixed K m to the value obtained from the literature on the assumption that drug affinity would not change much in the perfused rat liver. Then V max , K 12_r and K 21_r were optimized based on plasma concentration data and V max B, K m B and CL met were optimized by curve fitting the accumulated biliary excretion data [15]. Since the oral bioavailability of GL is quite low (,4% in rat and cannot be detected in human) [5,16,23,24], the absorption of GL after oral administration was ignored. The hydrolysis rate constant of GL (KH) in the large intestine was taken as the value determined by incubation of GL in gastrointestinal contents in vitro [23]. Other physiological parameters included in the gastrointestinal model were taken from the literature and are listed in Table 2 [17].

Description of the GA model
As in the GL model, four tissue compartments and one blood compartment were modeled to describe GA disposition in rat (shown in Figure 2A). The kidney and gut compartments were as described for the GL model. As nearly 100% of the radioactivity in a dose of tritium-labeled GA was found to be excreted into bile (mainly as metabolites) [25], it is clear that the liver plays an important role in the disposition of GA. GA can competitively inhibit the hepatic uptake of GL [26] indicating that the two compounds share the same sinusoidal transport system. When transported into the liver, GA is metabolized mainly by UDPglucuronosyltransferase and sulfotransferases in what are supposedly first order processes. The resulting phase II metabolites (GAM) are then excreted with bile into the intestinal tract and hydrolyzed to GA which then undergoes enterohepatic recycling.

Calibration of the GA model
On the assumption that GA shares the same hepatic uptake transporter as GL, the V max value which was considered to be related to the expression of the transporter was taken to be the same as that in the GL model leaving K m to be optimized. The metabolic clearance of GA by the human liver (CL in vivo, human ) was estimated based on its metabolism by human liver microsomes (HLM) in vitro [27] and further adjusted by the unbound fraction of drug in the incubation medium f u_inc using equations 2 and 3 where logD is the measure of lipophilicity of the compound [28].
The metabolic clearance of GA by the rat liver (CL int vivo, rat ) was then estimated by scaling on the basis if liverweight as given in equation 4: The absorption rate constants of GA in each intestinal segment were obtained by the in situ loop method provided by Ploeger et al [8,24]. The efflux clearance of GA from the liver to the venous blood (CL eff ) was assumed to be 15.6-fold greater than that of GL based on the fact that the passive permeability of GA in the gut lumen was estimated to be 15.6-fold greater than that of GL [24]. The data set of concentration-time profiles after i.v. administration of different doses of GA to rats with bile fistulas was used to optimize the parameters K m , K 12_r and K 21_r [29]. Like GL, the biliary excretion of GAM is via Mrp2 [5] for which V max B and K m B were assumed to be the same as those of GL and was validated using experimental data for: (1) The accumulated biliary excretion of 3MGA after i.v. administration of 3MGA at 5 mg/kg [5]; and (2) the accumulated biliary excretion of GAM after i.p. administration of GA at 25 mg/kg [25]. Then the colonic hydrolysis rate constant of GAM (KH GAM ) was optimized according to the plasma concentration-time profiles of GA after i.v. administration of GA at 5.7 mg/kg in rats without bile fistulas [30]. All other parameters were obtained from the literature and are listed in Table 3.
Equations for the PBPK models of rat and human are listed in detail in the Appendix S1.

Development of human PBPK model for GL and GA: from rat model to human model
The structure of the model of GL and its metabolites in human is similar to that in rat shown in Figure 2B. As the absorptive properties of the compounds in the different regions of the gastrointestinal tract are unclear, we divided the gastrointestinal tract into only three parts viz the stomach, small intestine and large intestine. The physiological parameters in human are listed in Table 2. Moreover, in the human model, GL and GAM are stored in the gallbladder from which they undergo biliary excretion into the intestinal lumen after the ingestion of fat at meal times [9]. The rate constants for hydrolysis of GL and GAM in the gut lumen were calculated based on the literature [32] and K 12_r and CL met in the human model were scaled from those in the rat model based on the cardiac output and liver weight respectively. Tissue to plasma partition coefficients for kidney (P k ) and gut (P g ) in human were assumed to be the same as those in rat.
The initial values of CL up and CL b in the human model were extrapolated from the corresponding clearances in the rat model using the allometric equations 5 and 6 where BW rat and BW human are the body weights of rat (250 g) and human (70 kg) respectively: The absorption rate constant of GL during enterohepatic recycling in human was obtained from the literature [19]. CL up and CL b of GL were then optimized according to plasma concentration-time data for a single 120 mg i.v. dose of GL. The simulated profiles following other i.v. doses of 40 and 80 mg [16] were compared with the observed data to further validate the model. CL up and the absorption rate constant (K abs ) of GA were optimized according to plasma concentration-time data for a single 130 mg p.o. dose [10]. The values were further validated for the following scenarios: Multiple p.o. administration of GA at 130 mg/day for 5 days [10]; a single p.o. administration of 225 mg GL; and a single p.o. administration of licorice equivalent to 225 mg GL [9]. For the licorice dose, the gut microbiota-mediated hydrolysis rate constant (KH) was assumed to be 4-fold greater since other components in licorice increase the extent of hydrolysis as shown in an in vitro study in rat [23].

Development of the pharmacodynamic model in human
As mentioned previously, GA reversibly inhibits 11b-HSD 2 in the kidney and gives rise to pseudoaldosteronism. The scheme of the PD model for the 11b-HSD 2 associated renin-angiotensin-aldosterone electrolyte system is shown in Figure 2C. The 11b-HSD 2 module was first developed from the Ploeger model with some modification. The reversible inhibition of the conversion of cortisol into cortisone by GA can be described by the previously developed inhibitory Imax model [10], where C is the concentration of GA in the kidney and IC50 is the concentration producing 50% of the maximum inhibition. The initial value of IC50 was assumed to be equal to the corresponding value found in rat kidney microsomes i.e. 0.32 mM [33]. Kidney exposure to cortisol was described by equation 8: where F k is the kidney concentration of cortisol in human, k 0 is the rate of cortisol formation (assumed to be constant), V k is the kidney volume and k ox and k ur_F are the 11b-HSD 2 metabolic clearance and urinary clearance of cortisol, respectively. For simplicity, we assumed that the elimination of cortisol takes place by these two pathways only. We also assumed that the change in cortisone concentration in plasma was only due to its formation by 11b-HSD 2-mediated oxidation of cortisol and urinary excretion. Accordingly, the change in plasma concentration of cortisone was described by equation 9: where k ur_E is the urinary clearance and V app is the apparent distribution volume according to the plasma cortisone concentration. As the urinary cortisol:cortisone ratio is the generally accepted biomarker of 11b-HSD 2 inhibition, the urinary content of cortisol and cortisone were calculated as follows: and the 24 h urinary cortisol:cortisone ratio (R) was described as: The parameters k ur_F and k ur_E were determined from the baseline 24 h urinary excretions of cortisol and cortisone [34] using the following equations: k ur E (L=h)~u rinary excretion rate(nmol=day)=24 where F k0~F0 NK tp cortisol ð15Þ F 0 and E 0 are the baseline plasma concentrations of cortisol and cortisone respectively, F k0 is the baseline kidney concentration of cortisol and K tp_cortisol is the kidney:plasma concentration ratio of cortisol [35]. When GL is absent and equilibrium has been reached, equation 16 can be derived from equation 9: and equation 17 from equation 8: The IC50 and V app value were optimized based on the 24 h cortisol:cortisone ratio for 10 days after the first day of daily consumption of GA at 130 mg/day for five days.
Once 11b-HSD 2 is inhibited, the increased concentration of cortisol activates the MR and further induces sodium and water retention, hypokalemia, and hypertension. Ikeda et.al [36] developed a comprehensive network for overall regulation of body fluids which we used with some modification in the present study as the electrolyte module to predict the serum potassium and sodium levels. Moreover, several studies have reported that a significant suppression of the renin-angiotensin-aldosterone axis occurs after over-consumption of licorice and GA taken to be a compensatory response to the excessive activation of the MR by cortisol [34,37]. Guillaud and Hannaert [38] developed a computational model of the renin-angiotensin system which can output the plasma renin, angiotensin I and angiotensin II levels but not aldosterone level, and successfully predicted the effect of renin inhibitor aliskiren. We then developed the mathematical relationship between aldosterone and angiotensin II given in equation 18 and further linked it with the Ikeda model shown in Figure 2C. Since angiotensin II and potassium are primary secretagogues modifying aldosterone release [39], we assumed that generation and secretion of aldosterone are mainly regulated by the two factors shown in equation 18 (i.e. the generation of Aldo stimulated by Angiotensin II and potassium in Figure 2C). Here K gen_aldo and K deg_aldo are the zero-order release rate constant and first-order elimination rate constant of aldosterone, S max is the maximum stimulation, SC50 is the concentration producing 50% of maximum stimulation and r is the Hill coefficient. These parameters were obtained by curve fitting the experimental data [40,41]. The RAAS is further regulated by the arterial pressure, venous pressure and glomerular filtration rate (GFR) which can be outputs of the Ikeda model [38].
The body fluids model developed by Ikeda consists of 7 blocks, 5 of which were the main focus of the present study and are described briefly as follows: Block 1, the cardiovascular system, calibrates the relationship between blood volume and systemic arterial and venous pressure; Block 3, the extracellular space, calibrates the changes in water in the vascular and extracellular space; Block 4, the intracellular space and electrolytes, calibrates the potassium and sodium concentrations in extracellular fluids; Block 6, the kidney, calibrates the urinary excretion of water, sodium and potassium; and Block 7, the controller of renal function, calibrates GFR. Block 2 and 5 correspond to the predictions of the respiratory system and urinary anion excretion respectively which are not the main outputs in our model. The Ikeda model was described mathematically as a set of nonlinear differential and algebraic equations of more than 200 variables. Behavior of the model for various kinds of inputs simulated with a digital computer was in good agreement with a number of experimental results pertaining to body fluids and electrolytes [36]. The scheme of the reninangiotensin-aldosterone-electrolyte network model modified on the basis of the original literature [36,38] is shown in Figure 2C with equations given in detail in Appendix S2. The effect of cortisol on the MR was added to the aldosterone effect in the electrolyte module as indicated in equation 19 (i.e. the stimulation of MR by Aldo and cortisol in Figure 2C): where ALD 0 is a normalized value presenting the effect of MR defined by Ikeda et al. and [Aldo] 0 and [F k ] 0 are the baseline concentrations of aldosterone in plasma and cortisol in kidney respectively. ALD 0 then affects the urinary excretion of potassium and sodium (YKU and YNU) in Block 6 of the Ikeda model. Finally the electrolyte module outputs the serum potassium and sodium levels.

Simulation software
All simulations were performed using MATLAB (the MathWorks Inc., Natick, MA, USA). The PKPD model was constructed as a set of ordinary differential equation (ODEs), the integration of which was performed using the fourth order Runge-Kutta method. The accuracy of the prediction was graphically evaluated by superimposing the concentration-time profile observed in vivo to the simulated ones.

PBPK modeling of GL and its metabolites in rat
In this study, a PBPK model was developed to predict the PK profiles of GL and GA after i.v. and p.o. administrations. Drug specific parameters used to develop the semi-PBPK model are listed in Table 3. Optimization of the parameters K 12 , K 21 , V max , V max B, K m B and CL met in the GL model according to plasma concentration-time profiles and accumulated biliary excretion are shown in Figure 3A-B. It can be seen that the predicted lines give a good fit to the observed data. The optimized value of V max (8.96 mmol/h) is quite close to that (10.1 mmol/ h) calculated from isolated liver perfusion experiments, indicating that fixing K m to the experimental value was a reasonable strategy. We further simulated different scenarios of i.v. administration of GL [18,19,31] and validated the model by comparing the simulations with data reported in the literature ( Figure 3C-D). The results show that the simulated and reported profiles are comparable.
Curve fitting of K m , k 12 and k 21 in the GA model to concentration-time profiles for i.v. administration of 2-20 mg/kg GA to rats with bile fistulas is shown in Figure 4A. Curve fitting of the 3MGA plasma concentration-time profile and subsequent simulation of its biliary excretion after i.v. administration of 3MGA at 5 mg/kg [5] and of GAM biliary excretion after i.p. administration of GA at  [15], (C) 10 mg/kg [13], (D) 25 and 100 mg/kg [18,19,31]; accumulated biliary excretion of GL at (B) 5-50 mg/kg [15] and (C) 10 mg/kg [13]. Experimental data (fitset and testset) are shown as symbols; the lines represent the prediction of the GL model. 25 mg/kg [25] are shown in Figure 5. After optimization of KH GAM by curve fitting of plasma concentration-time data for GA administration (i.v. 5.7 mg/kg [30]), the GA model was further validated for the following scenarios: 60 mg/kg i.v. [18] and 5.7 mg/kg p.o. GA; 10 mg/kg p.o. [30] and 100 mg/kg p.o. GL (experimental studies) shown in Figure 4B-C. Finally we simulated the kidney and plasma exposure to GA after p.o. administration of 200 mg/kg GL (experimental study) and compared the results with the observed data as shown in Figure 4D. The results show that the model predicts well the kidney exposure to GA which guarantees a reliable input into the PD model.

PBPK modeling of GL and its metabolites in human
The estimated PK profiles following three i.v. doses of GL were comparable with literature data as shown in Figure 6A. In the GA model, the CL up value and rate constants for absorption in the small intestine and colon were optimized based on the plasma concentration-time profile for p.o. administration of GA ( Figure 6B). Simulation of the PK profiles of GA for the three scenarios of administration mentioned in the Methods ( Figure 6B-C) illustrates that the model gives good predictions for the different scenarios. The drug specific parameters of the human PBPK model are listed in Table 4.

PK/PD modeling for prediction of the adverse effects of GL in human
The human PK/PD model was constructed in two stages: the first to predict the urinary cortisol:cortisone ratio (11b-HSD 2 module), the second to estimate the levels of serum potassium and sodium, plasma aldosterone etc. in the RAAS and electrolyte modules. Figure 7A shows the results of optimization of the IC50 and V app values based on the observed 24 h urinary cortisol:cortisone ratio after p.o. administration of 130 mg/day GA for five days. The optimized parameters are listed in Table 5. Validation and simulation of the urinary cortisol:cortisone ratio were further carried out for the following two scenarios: (1) Multiple p.o. administrations of 250 mg GA twice a day for 7 days ( Figure 7D); and (2) multiple p.o. administrations of 170 mg GA three times a day for 2 days (Table 6). In the first scenario, we also simulated the 24 h urinary excretion of cortisol and cortisone and compared the results with the observed data ( Figure 7B-C). It was found that the predicted and observed data were comparable indicating the resulting model was robust. Curve fitting the aldosterone model with the experimental data is shown in Figure 8A-B and the optimized parameters are listed in Table 5.
We then predicted the changes in the RAAS and serum electrolyte levels after p.o. administration of GA and licorice. The results shown in Tables 7 and 8 suggest that the RAAS was significantly inhibited presumably as a compensatory Figure 5. Plamsa concentration and accumulated biliary excretion of GAM. Left, plasma concentration and accumulated biliary excretion of 3MGA after i.v. 5 mg/kg 3MGA in rats [5]; Right, accumulated biliary excretion of GAM after i.p. 25 mg/kg GA in rats [25]. Experimental data are shown as symbols; the lines represent the predictions of the GA model.

Identification of the important factors relating to GL-induced pseudoaldosteronism
We further investigated the effects of changes in sinusoidal transport function, colonic transit time and activity of 11b-HSD 2, all of which have been reported to display high inter-individual variability, on the PK of GA and serum potassium level (chosen as the indicator of pseudoaldosteronism) after p.o. administration of 200 mg/kg GL in human.
As mentioned above, hepatic clearance is the major elimination pathway of GA and, according to the parameters in the PBPK model of GA, hepatic sinusoidal uptake clearance (20000 l/h) is almost 1000 times greater than sinusoidal efflux clearance (28.18 l/h) meaning sinusoidal uptake rather than intrinsic hepatic metabolism is the determinant of hepatic disposition. Moreover, expression of the hepatic sinusoidal OATP transporters shows a large inter-individual variability  suggesting that hepatic sinusoidal uptake function may contribute to the interindividual variability in GL-induced adverse effects. CL up is proportional to the V max of the uptake process which is further supposedly proportional to the expression of the sinusoidal transporter. OATP1B1 and OATP1B3 are the major transporters mediating GL and GA sinusoidal uptake in human and regulation of their expression, according to previous reports, is coordinated within an individual [42]. Thus, in the present study, the change in expression of OATP1B1 was taken as a measure of the change in CL up . Details are as follows. According to a previous report [43], there is a 4.9-fold inter-individual variability in OATP1B1 mRNA expression level in the liver. Thus the median was set to the normal value N and two other values of CL up (N/2.21 and 2.21*N) were chosen. The results ( Figure 9A-B) show that a decrease in CL up can significantly increase exposure to GA and subsequently decrease the serum potassium level, indicating that a decrease in sinusoidal transport function increases sensitivity to GL-induced adverse effects.
Since the adverse effect of GL is especially severe in elderly patients [44,45] and since, according to previous reports, age is an important factor impacting colonic transit time and the activity of 11b-HSD 2 [46,47], the effects of colonic transit time and activity of 11b-HSD 2 on the serum potassium level were further investigated. According to the reported values for colonic transit time and cortisol:cortisone ratio in the elderly, we set K co and k ox,0 to 0.0217 h 21 and 0.02 l/h respectively (the normal value of k ox,0 is about 1.8 times greater), predicted the plasma exposure to GA and potassium levels for one week and further compared them with those of using normal values of K co and k ox,0 . The results shown in Figure 9C-E suggest that a prolonged colonic transit time and a decrease in 11b-HSD 2 activity make subjects more susceptible to experiencing adverse effects and may explain why GL-induced pseudoaldosteronism mainly occurs in the elderly.
According to the literature mentioned above [43,46,47], the distribution of CL up in the general population and K co and k ox,0 in the elderly can be estimated and are listed in Table 9. Using the Monte Carlo method, we further simulated the dose limit for 1000 virtual elderly subjects according to the known distribution of the above three parameters in the elderly. Based on the normal range of serum potassium of 3.5-5.0 mmol/l, the median value (4.18 mmol/l) was chosen as the typical normal value of serum potassium. A series of doses were simulated for each subject and the individual dose limit leading to hypokalemia was determined as the dose giving a serum potassium level of 3.5 mmol/l on the 28 th day. The results are shown in Figure 10A. The distribution was well described by a logarithmic normal distribution and the distribution parameters m and s were found to be 5.6348 and 0.5449 respectively. Moreover, in the left part of the cumulative probability curve shown in Figure 10B, there is a critical value point of around 100 mg representing the highest rate of change in the cumulative probability. We further calculated it as the maximum value in the second derivative of the cumulative probability distribution function and obtained a value of 101 mg ( Figure 10C). This indicates that the risk of adverse effects increases abruptly above 101 mg making this the critical dose limit of GL causing hypokalemia in the elderly with a probability of 3.07%. We further estimated the probability of hypokalemia in the elderly for the upper dose limit recommended by SCF (100 mg/day) as 3% and DNIB (200 mg/day) as 27%. The results show that, although the dose is only doubled, the probability is unexpectedly 9 times higher. It is suggested that the risk of the adverse effect at the DNIB recommended dose (200 mg/kg) is substantially greater probably due to high popularity of licorice containing foods in the Netherlands.

Discussion
Consumption of GL in large amounts or for a long period of time leads to the adverse effects of pseudoaldosteronism which can be quite severe if ignored. In the present study, a semi-PBPKPD model was developed to predict the adverse effects of GL and facilitate its safe use both clinically and as a food additive.
Pseudoaldosteronism is caused by inhibition of 11b-HSD 2 by GA, the major metabolite of GL. Since the kidney is the target tissue for the adverse effects of GL, kidney exposure to GA is likely to be more important than systemic exposure and could be directly linked with the PD model. We made this prediction and validated it in the rat model and assumed that the plasma to kidney distribution in rat and human was the same. As a consequence, a semi-PBPK model was first developed in rat mainly to validate the model structure, clarify the mechanisms of action and provide reasonable parameters for the human model. The present liver model was developed mainly on the basis of Ishida's isolated perfused liver model in which the parameters of GL hepatic disposition were calculated from experimental data and used as initial values in optimization. The structure of the PBPK model can be simplified or complicated according to the objective based on a reasonable mechanism. According to our objective of identifying a high risk population in which sinusoidal transport may be an important factor causing variability in adverse effects, the present model connected sinusoidal uptake, intrinsic hepatic metabolism and biliary excretion in a tandem pattern allowing analysis of the impact of sinusoidal transport on exposure to GL and GA. In contrast, in Ploeger's PBPK model, partial sinusoidal uptake and biliary excretion were combined into one process, making it difficult to further analyze the rate sensitive process in drug elimination.
Simple exposure-response PD models cannot explain the complicated mechanisms involved in such things as signal transduction and feedback loops making it difficult to predict clinical indicators far from the site of action. A mechanistic model can, however, include factors and processes involved in wellclarified mechanisms of action making it possible to simulate clinical reality, to predict more than one indicator at the same time and to establish the quantitative relationship between target-related and end point-related indicators. Thus, in order to achieve this objective, we developed a mechanism-based PD model. The PD model was further divided into three parts: the 11b-HSD 2 module (target The observed value [37] is expressed as mean (SD).
doi:10.1371/journal.pone.0114049.t008  related), the RAAS module (feedback related) and the electrolyte module (end point related). The 24 h urinary cortisol:cortisone ratio is a generally accepted biomarker of the inhibiton of 11b-HSD 2 and sensitive to changes in the GA level in the kidney. In the present PD model, not only the ratio but also the urinary excretion of cortisol and cortisone could be predicted and further validated by comparison with experimental data. This made the further prediction of kidney exposure to cortisol more reasonable. In addition to the 24 h urinary cortisol:cortisone ratio, serum electrolytes (potassium and sodium) and RAAS are also important indicators of GL-induced pseudoaldosteronism and are more indicative of the occurrence of adverse effects. In the present PD model, the plasma levels of angiotensin II, aldosterone, potassium and sodium were estimated and compared with the normal ranges of each indicator to predict the occurrence of adverse effects. Since there is no gold standard for GL-induced pseudoaldosteronism and since, of the various indicators of pseudoaldosteronism, serum potassium (hypokalemia) is an independent risk factor with the narrowest normal range (3.5-5 mmol/l) [48] and quite sensitive to the over-activation of MR [34,37], we chose serum potassium as the main indicator in the simulation and prediction of the safe dose of GL. After validation of the PKPD model, it was used to analyze the safe dose of GL. According to the SCF, a prudent dose of GL not exceeding 100 mg/day is recommended which is the LOAEL value in the most sensitive individuals. The simulation of the population mean in this scenario by the present model shows that 100 mg/day produces a less than 0.1 mmol/l decrease in the serum potassium level, which is considered safe for normal people. The screening results show that sinusoidal transport function, colonic transit time and the activity of 11b-HSD 2 are all important factors contributing to the inter-individual variability in experiencing adverse effects. In the present study, we considered the elderly as a high risk population and estimated the distribution of the individual dose limit based on data from recent studies.
In reality, describing an adverse effect should include three major factors: the dose limit, the specified population and the corresponding incidence of adverse effects. The present study provided the probability distribution of the individual dose limit in a virtual elderly population and calculated the critical value as the dose limit recommended by the present model. Interestingly, the critical value of 101 mg in our virtual elderly population was almost identical with the dose limit of 100 mg recommended by SCF. However, it must be pointed out that the SCF dose limit is an experimental value based on the clinically observed symptoms of pseudoaldosteronism while our critical value is based only on the serum potassium level. Also, the PKPD model provides the probability of increasing dose in a virtual population of 3.07% at the critical value which cannot be obtained from clinical observations of adverse effects due to ethical concerns. This is the most important advantage of modeling and simulation compared to clinical experiments. In further clinical trials, if new factors such as drug-drug interactions or the effects of polymorphism, or other metabolites of GL are discovered to be involved in producing adverse effects, they can be introduced into the model. In addition, the model can be used to simulate scenarios which may not be included in clinical trials.
In summary, we have developed a comprehensive PBPK/PD model based on PK data and pharmacology to predict the occurrence of GL-induced pseudoaldosteronism. In this model, not only the urinary cortisol:cortisone ratio but also RAAS and serum potassium and sodium levels could be predicted. Furthermore, the present model simulated the probability distribution of the individual dose limit in a virtual population and provided the critical value of the dose limit to be recommended when GL is used as a food additive. This work would be useful to improve the safe use of GL as a food additive and drug.

Supporting Information
Appendix S1. Main equations of the PBPK model for GL and its metabolites in rat and human.