Mathematical Models in the Description of Pregnane X Receptor (PXR)-Regulated Cytochrome P450 Enzyme Induction

The pregnane X receptor (PXR) is a drug/xenobiotic-activated transcription factor of crucial importance for major cytochrome P450 xenobiotic-metabolizing enzymes (CYP) expression and regulation in the liver and the intestine. One of the major target genes regulated by PXR is the cytochrome P450 enzyme (CYP3A4), which is the most important human drug-metabolizing enzyme. In addition, PXR is supposed to be involved both in basal and/or inducible expression of many other CYPs, such as CYP2B6, CYP2C8, 2C9 and 2C19, CYP3A5, CYP3A7, and CYP2A6. Interestingly, the dynamics of PXR-mediated target genes regulation has not been systematically studied and we have only a few mechanistic mathematical and biologically based models describing gene expression dynamics after PXR activation in cellular models. Furthermore, few indirect mathematical PKPD models for prediction of CYP3A metabolic activity in vivo have been built based on compartmental models with respect to drug–drug interactions or hormonal crosstalk. Importantly, several negative feedback loops have been described in PXR regulation. Although current mathematical models propose these adaptive mechanisms, a comprehensive mathematical model based on sufficient experimental data is still missing. In the current review, we summarize and compare these models and address some issues that should be considered for the improvement of PXR-mediated gene regulation modelling as well as for our better understanding of the quantitative and spatial dynamics of CYPs expression.


Introduction
Nuclear receptors (NRs), which are ligand-or hormon-activated transcription factors, are essential in the organism in many processes, such as hormonal regulation, homeostasis of endogenous compounds and ions, the regulation of both intermediary and xenobiotic metabolism, and the regulation of cell differentiation, development, and proliferation [1]. Furthermore, these receptors are associated with numerous health conditions and unbalanced physiological processes, for instance cardiovascular diseases, inflammation, neurodegenerative or liver deseases, and cancer. Highly relevant diseases and In the current review, we summarize quantitative models that describe PXR-mediated regulation of cytochrome P450 3A gene expression with several experimental models for primary human or rat hepatocyte cell lines. In particular, we focus on mechanistically and biologically based cellular models that concentrate on molecular aspects of PXR-mediated regulation of the most important target gene, CYP3A4, or its rodent orthologues.

Pregnane X Receptor Characterization
The human PXR is the product of the NR1I2 gene, which is located on chromosome 3, locus 3q11-q13. 3. The NR1I2 gene comprises 10 exons separated by nine intronic regions [22,23]. Several protein isoforms of the human PXR have been characterized with differential transcriptional activity. The transcript variants 1 (PXR1) and 2 (PXR2) respond to an agonist by activating target gene expression, while transcript variants 3 (PXR3) and the truncated variant 4 (PXR4) do not induce target gene expression [24]. The identification of PXR as a xenosensor presents a molecular basis for the species specificity of the induction of the CYP3A genes [25]. High homology (95% at the amino acid level) has been reported between human PXR (hPXR) and mouse PXR (mPXR) in the DNA-binding domain (DBD) [26]. Thus, they can share PXR binding sites found in promoters of both the human CYP3A (CYP3A4, CYP3A5, and CYP3A7) and rodent Cyp3a (e.g., Cyp3a11) genes. Nevertheless, the homology is much lower (73% at the amino acid level) in the ligand-binding domain (LBD), which may illustrate the ligand specificity between these receptors. Furthermore, the X-ray crystal structure analysis of the PXR LBD supports this notion [27]. There are significant variations between PXR species regarding PXR activators. The ligand specificity of hPXR is very broad, mainly for large ligands, such as rifampicin, whereas mPXR has a narrower ligand specificity. For instance, mPXR is more sensitive to pregnenolone 16α-carbonitrile than hPXR, while hyperforin and rifampicin can activate hPXR more than mPXR [28]. Additionally, it has also been reported that dexamethasone is a more potent ligand for mPXR than for the human receptor [29,30]. PXR has a spherical ligand-binding pocket, which has been shown to be at least twice as large as those of the other steroid hormone or retinoid receptors [27]. The PXR ligand-binding pocket was found to be hydrophobic and flexible. Therefore, these structural properties may have accounted for the high ability of this receptor to recognize a diverse range of xenobiotics [27]. The species origin of the PXR receptor, rather than the promoter structure of CYP3A genes, dictates the species-specific pattern of CYP3A inducibility [31]. These results have encouraged the creation of "humanized" hPXR transgenic mice, where the human counterpart hPXR genetically replaces the mouse PXR in the liver [26].

PXR Localization
PXR expression occurs mainly in the liver, intestine, and to a lesser degree in the kidney. Minor expression of hPXR/mPXR mRNAs may appear in other tissues, such as the lung, stomach, peripheral blood monocytes, uterus, ovary, breast, adrenal gland, bone marrow, and some regions of the brain [11,16].
The mPXR is primarily located in the cytosol of untreated liver cells, where mPXR binds with the cytoplasmic CAR retention protein (CCRP) and the heat shock protein 90 (Hsp90) to form a protein complex [11]. The latter retains the cytosolic localization of PXR. When the ligand binds to mPXR, the PXR dissociates from the multi-protein complex and translocates to the nucleus in mouse hepatocytes to regulate gene transcription [32,33] (Figure 1). Nevertheless, nuclear localization of hPXR has been found in mammalian tumor-derived cell lines [34].

PXR Transcriptional Machinery
As a typical NR, PXR has the DBD at the N-terminus and the LDB at the C-terminus. The DBD is capable of binding to regulatory DNA sequences [14]. Typically, there are two transcriptional activation domains in a nuclear receptor: the activation function 1 (AF-1), which resides in the N-terminal domain, and the AF-2, which is present in the C-terminal portion of the LBD. In the case of PXR, only one AF-2 is present in the LBD domain. The LDB has many functions, including ligand binding, dimerization, transcriptional activation, and interaction with transcriptional cofactors. AF-2 is the helix at the C-terminus that controls transcriptional activation through coactivators by conformational rearrangements of LBD or gene suppression through transcriptional corepressors recruitment [35,36].
The activation of nuclear receptors mostly takes place in the cytoplasm, where unliganded nuclear receptors typically reside in the complex with chaperones (such as Hsp90) and with co-chaperone CCRP. The schematic of PXR nuclear receptor activation is shown in Figure 1. Xenobiotics bind to the LBD of the PXR after they enter the cell. The binding of a ligand to the LBD results in a conformational change in the AF-2 that disrupts interactions with transcriptional corepressor proteins, such as NCoR and SMRT, triggers the release of the co-chaperone complex, and permits cytosolic-nuclear translocation of liganded PXR along cytoskeletal tracks toward the nucleus. Since the molecular weight of PXR is about 50 kDa, the nucleocytoplasmic transport is specifically regulated by an energy-dependent reaction (Nigg, 1997;Mattaj and Engelmeier, 1998). In the nucleus, formation of heterodimers with another NR called Retinoid X receptor-α (RXRα) and interactions with transcriptional coactivator proteins, such as members of the p160/steroid receptor coactivator (SRC) family, take place [8,33,37]. Afterwards, the PXR/RXRα heterodimer binds to the DNA through a pregnane X receptor response element (PXRE) and promotes the transcription of the target genes' mRNAs [15,38]. These mRNAs then move to the cytoplasm, where they are subjected to translation into proteins which, as a result, do their function in metabolizing xenobiotics that invaded the cell [39].
In the review by Moore et al. [7], the authors reviewed some of the major coactivators and corepressors involved in the PXR regulation of drug-metabolizing enzymes and transporters. Steroid receptor coactivators 1 (SRC1/NCOA1) and 2 (SRC2/GRIP1), the NR interacting protein 1 (NRIP1/RIP140), the peroxisome proliferator-activated receptor gamma coactivator α (PGC-1α), and the Forkhead transcription factor FKHR (FOXO1) were identified as PXR coactivators [40]. Small heterodimer partner (SHP/NC0B2) and NR corepressor 2 (NCoR2/SMRT) were determined to be corepressors [40][41][42]. According to [41], in HepG2 cells, SHP hinders PXR binding to DNA in a ligand-dependent manner. Signaling cascades play a role in regulating PXR activity. In particular, the stimulation of PXR-mediated CYP3A induction by protein kinase A (PKA) has been shown [43] and PXR was reported to be phosphorylated by PKA in vitro. However, protein kinase C (PKC) was found to be a PXR signaling repressor [8,44]. The degradation of PXR via proteasomes is targeted by conjugation with a polyubiquitin chain. This "ubiquitination" process is governed by a sequential pathway of three different enzymes (E1, E2, and E3) and results in covalent binding of ubiquitin to PXR [45]. Other papers reported the involvement of acetylation in regulating PXR function [46,47]. Acetylation represents one of the major mechanisms for post-translational modification of proteins, including PXR [48]. PXR can also be regulated by other different signaling pathways, such as SUMOylation [8]. In the regulation, two feedback loops are involved. First, PXR triggers CYP3A4 mRNA transcription resulting in CYP3A4 protein synthesis in the cytosol. Subsequently, cytosolic CYP3A4 metabolically inactivates rifampicin, a substrate of CYP3A4. In the second feedback loop, activated PXR downregulates its own expression via the PXR (NR1I2) gene promoter. A potential third proposed feedback loop (not shown) supposes that PXR downregulates Small heterodimer partner (SHP) expression. SHP is a corepressor of PXR [49]. (1) passive diffusion of a ligand into cells; (2) binding to PXR that forms a cytoplasmic complex with chaperones (Hsp90), CCRP, and SMRT. The binding to a ligand promotes release of these proteins from PXR; (3) cytoplasm-nuclear translocation; (4) binding of coactivators and other nuclear receptors (such as RXRα) to form a transcriptional complex; (5) binding of the complex to PXR response elements in target genes' promoter regions; (6) recruitment of transcription machinery factors and mRNA synthesis; (7) export of mRNA into cytosol and translation. HNF4α, hepatocyte nuclear factor 4α, a nuclear receptor that coactivates human PXR (hPXR); PGC1α, Peroxisome proliferator-activated receptor gamma coactivator 1-alpha, a coactivator of PXR; PolII, RNA polymerase II; TBP, TATA-binding protein; TFII, transcription factor II D.

PXR Target Genes Regulation
PXR is currently regarded as a master transcription factor for xenobiotic-and drug-inducible expression of key genes that encode members of the phase I and phase II metabolic enzymes and drug transporters. Additionally, previous reports suggest that PXR plays an integral role in endobiotic metabolism by regulating important genes which have major roles in glucose, lipid, and bile acid metabolism [11]. A recent paper documented that PXR upregulates 164 genes, but downregulates the expression of 334 genes, in primary human hepatocytes [50].

CYP3A4 Gene Regulation via PXR
Trans-as well as cis-regulator transcriptional elements involved in CYP3A4 gene expression have been comprehensively studied in the past. Apart from the basal promoter with an ER6 response element, the distal xenobiotic responsive enhancer module (XREM) is the key gene expression enhancer controlled by PXR. Depending on the species, these elements consist of DR3 motifs or ER6 motifs in the CYP3A4 gene promoter and function in a tissue-specific manner [15,51].

PXR-Mediated Regulation of Other Target Genes
In addition to the CYP3A4 gene, PXR was characterized as an important metabolic regulator for the induction of the CYP2B6, CYP2C8/9, CYP2C19, CYP3A4/5/7, and CYP2A6 cytochrome P450 enzymes, phase II enzymes, such as UDP-glucuronyltransferases and sulfotransferases, and the upregulation of transporters, such as ABCB1/P-glycoprotein.

Figure 1.
Schematic picture of pregnane X receptor (PXR)-mediated regulation of CYP3A4 gene expression. In the regulation, two feedback loops are involved. First, PXR triggers CYP3A4 mRNA transcription resulting in CYP3A4 protein synthesis in the cytosol. Subsequently, cytosolic CYP3A4 metabolically inactivates rifampicin, a substrate of CYP3A4. In the second feedback loop, activated PXR downregulates its own expression via the PXR (NR1I2) gene promoter. A potential third proposed feedback loop (not shown) supposes that PXR downregulates Small heterodimer partner (SHP) expression. SHP is a corepressor of PXR [49]. (1) passive diffusion of a ligand into cells; (2) binding to PXR that forms a cytoplasmic complex with chaperones (Hsp90), CCRP, and SMRT. The binding to a ligand promotes release of these proteins from PXR; (3) cytoplasm-nuclear translocation; (4) binding of coactivators and other nuclear receptors (such as RXRα) to form a transcriptional complex; (5) binding of the complex to PXR response elements in target genes' promoter regions; (6) recruitment of transcription machinery factors and mRNA synthesis; (7) export of mRNA into cytosol and translation. HNF4α, hepatocyte nuclear factor 4α, a nuclear receptor that coactivates human PXR (hPXR); PGC1α, Peroxisome proliferator-activated receptor gamma coactivator 1-alpha, a coactivator of PXR; PolII, RNA polymerase II; TBP, TATA-binding protein; TFII, transcription factor II D.

PXR Target Genes Regulation
PXR is currently regarded as a master transcription factor for xenobiotic-and drug-inducible expression of key genes that encode members of the phase I and phase II metabolic enzymes and drug transporters. Additionally, previous reports suggest that PXR plays an integral role in endobiotic metabolism by regulating important genes which have major roles in glucose, lipid, and bile acid metabolism [11]. A recent paper documented that PXR upregulates 164 genes, but downregulates the expression of 334 genes, in primary human hepatocytes [50].

CYP3A4 Gene Regulation via PXR
Transas well as cis-regulator transcriptional elements involved in CYP3A4 gene expression have been comprehensively studied in the past. Apart from the basal promoter with an ER6 response element, the distal xenobiotic responsive enhancer module (XREM) is the key gene expression enhancer controlled by PXR. Depending on the species, these elements consist of DR3 motifs or ER6 motifs in the CYP3A4 gene promoter and function in a tissue-specific manner [15,51].

PXR-Mediated Regulation of Other Target Genes
In addition to the CYP3A4 gene, PXR was characterized as an important metabolic regulator for the induction of the CYP2B6, CYP2C8/9, CYP2C19, CYP3A4/5/7, and CYP2A6 cytochrome P450 enzymes, phase II enzymes, such as UDP-glucuronyltransferases and sulfotransferases, and the upregulation of transporters, such as ABCB1/P-glycoprotein.
In the paper by Goodwin [52], it was shown that PXR is capable of activating the phenobarbital-responsive enhancer module region of the CYP2B6 gene. CYP2C8 and CYP2C9 are regulated by PXR via common NR AGGTCA-based DNA response element motifs, including everted repeats spaced by six base pairs (ER6) [53]. The P-glycoprotein transporter seems to be regulated by PXR through a DR4 motif in the upstream enhancer at about −8 kilo base pairs from transcription start [54]. In addition to the P-glycoptotein transporter, multidrug resistance protein 2 (MRP2, ABCC2), multidrug resistance-related protein-3 (MRP3, ABCC3), and organic anion transporter polypeptide-2 (human SLC22A7, mice oatp2) are regulated by PXR activation [54][55][56][57].

Brief Overview of Physiological Functions of PXR
PXR is today known to have numerous functions for various physiological purposes. Regulation of detoxification enzymes also contributes to the elimination of toxic endogenous compounds. Bile acids are a family of endogenous PXR ligands identified shortly after PXR had been cloned. Accumulation of bile acids in the liver can cause hepatocyte toxicity, and impaired elimination of bile acids from the body can be accompanied with pruritus during cholestatic diseases. PXR has been reported to function as a lithocholic acid sensor and has a crucial role in the detoxification of cholestatic bile acids [62,63]. Therefore, rifampicin, a prototype PXR ligand, is now used to alleviate pruritus by augmentation of bile acids metabolism and elimination [64]. Accumulation of bilirubin in the blood may cause neurotoxicity. PXR has been shown to stimulate the expression of a variety of essential factors in the bilirubin clearance pathway [65,66].
In addition, PXR plays an important endobiotic role in vitamin K homeostasis [67] and in vitamin D catabolism through the upregulation of CYP3A4 and CYP24 [68,69]. Importantly, vitamin E [70,71], beta-carotene [72], and some other bile acids [62] have been proposed as putative endogenous PXR ligands, which indicates that the physiological function of the compounds could be partly mediated by PXR activation. Further work shows that activation of PXR by rifampicin decreases the activity of NF-κB, an essential regulator of inflammation and the immune response [73,74]. PXR has also been considered in tumor progression [75].
A broad area of physiological functions have been proposed for PXR with respect to lipid, cholesterol, and glucose metabolism regulation. For detailed information, refer to the recent reviews [11,76,77].

Compartmental Models
The main regulatory role of nuclear receptors, such as PXR, is up-and downregulation of gene expression resulting from augmented transactivation or trans-repression of its target genes. Gene expression is the outcome of the perturbation of a hierarchically organized and highly controlled network of interacting factors in signaling pathways and gene regulation networks in the cells, where ligand-activated nuclear receptors respond to an external stimulus. This process is distributed into several cellular compartments, including the nucleus and the cytosol. To understand these sophisticated networks, it is often beneficial to study the overall system by computational techniques based on complex biological knowledge, i.e., by a systems biology approach with a (semi-)mechanistic model rather than individual processes .
In a systems biology approach, the quantitative relationships in biologic systems are conveniently expressed by mathematical models. The important advantage of the models is that we can simulate different processes within different intervals, follow the dynamics of output parameters, and examine the effects of different parameters on different endpoints. Classical pharmacokinetic compartmental models, for example, describe the way one or several substances are distributed among distinct parts of a system through time. The system is conceptualized as consisting of one or more compartments, with the actual amount of each substance being represented in each compartment. The model then characterizes the transport of the substances from one compartment to another and/or the interaction of the substances inside a compartment, where any concentration inhomogeneity is ignored. Among the pharmacokinetic models, there are three main types: compartmental models, physiologically based models, and non-compartmental models [78]. Compartmental models represent the traditional and the most widely used approach in physiologically based pharmacokinetic (PBPK) and pharmacokinetic/pharmacodynamic (PK/PD) modelling [79].
The compartments are usually related to different sections of a body (such as organs) or different organelles of a cell type, each characterized by its volume of distribution (which may coincide with physiological volumes), within which the concentration of a drug is assumed to be homogeneous and to vary through time only. The mathematical description of the transport and interaction in time then consists of a system of ordinary differential equations (ODEs). Some of the parameters of the system of ODEs may be found in the literature and some donor-dependent parameters may be measured or derived directly. However, in practice, several parameters of the system of ODEs usually need to be estimated from in vitro or in vivo experiments through curve fitting. This can lead to several modelling stages: the parameters of a first simple model are fitted and the model is evaluated using obtained experimental data. The result may indicate where the model needs to be refined and an advanced model can be built with additional parameters, which in turn need to be fitted.
Importantly, the most useful technique for modelling the dynamics of gene network regulation over time employs compartmental modelling [80]. Among NRs, which regulate detoxification or metabolic enzymes, the glucocorticoid receptor (GR) seems to be the only receptor that has been extensively studied, and numerous mathematical models of several generations have been published to simulate its trans-activation in its target genes' regulation [81]. Notably, by employing mathematical models, both ligand-dependent downregulation and circadian oscillation of GR receptor expression can be simulated and predicted [82]. Very recently, a multiscale pharmacodynamic model was developed to characterize the overall GR-mediated effects on transcriptome, proteome, and enzymatic activities after a single dose of methylprednisolone [83]. In the work, the authors also considered post-transcriptional processes in their mechanistic mathematical model that control some GR target gene expression.
In addition to mechanistic and biologically-based models, indirect response models and PK/PD modelling have been used to describe the effect of ligands of some nuclear receptors. Recently, population PK-PD models have been used to describe the exposure-response relationship between plasma rifampicin and 4β-hydroxycholesterol, an endogenous product of CYP3A [84,85]. Similarly, an indirect PK/PD response model has been built for the peroxisome proliferator-activated receptor gamma (PPARγ) agonist pioglitazone in vivo [86]. In the NONMEM model, a plasma glucose lowering effect was selected as a surrogate for an anti-diabetic effect of the drug. In another clinical study, a PK/PD model has been built to model the side effects of LY500307, a highly selective estrogen receptor β (ERbeta) agonist [87].
In the case of PXR, only three quantitative mechanistic mathematical models of PXR-mediated human CYP3A4 gene regulation in the liver have been published in the literature that incorporate data from humans or from human models [88][89][90]. In addition, one further model described rodent PXR regulation of its target genes (CYP3A1/2) in rats or in primary rat hepatocytes with the rodent PXR ligands dexamethasone, lithocholic acid, and pregnenalone-16α-carbonitrile [91]. Finally, Kolodkin et al. recently published an in silico model of the dual effects of the stress hormone cortisol on both PXR and glucocorticoid receptors activation [89]. The following table summarizes the application characteristics of the published models and specifies the software used for in silico computations ( Table 2).
For the first model in the table, we display below the corresponding system of ordinary differential equations explicitly to give an impression of the arising type of equations. For the other models, the arising systems of ordinary differential equations are more complicated, with their explicit form sometimes given in the corresponding supplementary material or merely indicated with a diagram created by the used software. For the first model, the equations are: Here, the variables X ext (t), X int (t), PR(t), mRN A(t), and CYP3A4(t) denote the concentrations of, respectively, the xenobiotic outside the cell, the xenobiotic inside the cell, the PXR/RXR heterodimer, mRNA, and CYP3A4; and d(t) is the time-dependent dosing function. The parameters k imp , k exp , k dis , k mRN A , k mRN A,deg , k cyp , and k cyp,deg are first order constants for, respectively, the import and export of the xenobiotic, the heterodimer dissociation, the transcription rate of mRNA, the degradation coefficient of mRNA, the translation rate of CYP3A4, and the degradation coefficient of CYP3A4. k assoc is the association rate constant for the formation of PXR/RXR heterodimer, S PXR is the total system PXR concentration (bound and free), which is assumed to be constant, k met is a second-order metabolic constant, and p mRN A,back is the background production rate for mRNA. A more detailed description of the processes represented by the five equations, including an explanatory modelling figure, is given in the next section.
As seen above, the model considers PXR expression, which is denoted as S PXR . Importantly, this parameter can significantly vary in a human population due to the PXR expression variability. Highly variable expression of PXR mRNA and protein has been published many times before [92,93] and over 100-fold variability in PXR mRNA expression in both males and female can be found in larger human populations (GTEx Portal, Broad Institute of MIT and Harvard University). Numerous single-nucleotide polymorphisms of the NR1I2 gene have been described to affect PXR expression, having clinically significant manifestations or being associated with disease progression [94,95] (PharmGKB Database). In addition, different transcription variants can be present in different humans having different transcriptional activities [24]. Thus, the incorporation of the mechanistic and molecular aspects into a model can help us to reflect interindividual PXR expression and activity. The first quantitative (semi-)mechanistic biologically based compartmental model for PXR-mediated human CYP3A4 gene regulation has been described in the report by Luke et al. [90]. A schematic representation of the modelled processes is given in Figure 2, which can also be found in Luke et al. [90].

Model by Luke et al.
The first quantitative (semi-)mechanistic biologically based compartmental model for PXRmediated human CYP3A4 gene regulation has been described in the report by Luke et al. [90]. A schematic representation of the modelled processes is given in Figure 2, which can also be found in Luke et al. [90]. The model considers the standard processes connected with PXR signaling and PXR-mediated transactivation of its target genes, including the entry of rifampicin, a model PXR ligand, into hepatocytes by passive diffusion (described by Fick's law), ligand binding to PXR leading to the formation of PXR/RXR heterodimer in the nucleus, binding of heterodimer to DNA, mRNA transcription, and translation of mRNA to CYP3A4 proteins (described as a first-order process).
In addition, the model incorporates mRNA background production and degradation, CYP3A4 degradation, and, importantly, a feedback CYP3A4-mediated degradation of rifampicin. First-order kinetics are used to model the diffusion between compartments, heterodimer dissociation, mRNA transcription and translation, degradation of mRNA and CYP3A4, and the formation of the heterodimer, while the degradation of rifampicin by CYP3A4 is modelled by second-order kinetics. Incorporation of zero-order kinetics of mRNA background production and first-order kinetics of CYP3A4 mRNA degradation enables us to model steady-state CYP3A4 mRNA levels. This is very convenient, since there is a constitutive CYP3A4 mRNA as well as CYP3A4 protein expression in normal hepatocytes. In addition, this allows us to use fold-induction data to describe upregulation (induction) of CYP3A4 mRNA or protein from in vitro cellular experiments. These data can be easily obtained using standard qRT-PCR or Western blotting methods. In contrast, absolute quantification of mRNA transcripts is always difficult and needs sophisticated instrumentation or calibration.
The corresponding system of ordinary differential equations for the model is depicted above. Two forms of the model have been proposed in the work. The first one describes the situation in in vitro primary human hepatocytes in culture. It is composed of two compartments: the culture medium outside the plated cell culture (where rifampicin is present) and the interior of the plated cells. The model does not consider the compartmentalization of a cell into a nucleus and cytoplasm, The model considers the standard processes connected with PXR signaling and PXR-mediated transactivation of its target genes, including the entry of rifampicin, a model PXR ligand, into hepatocytes by passive diffusion (described by Fick's law), ligand binding to PXR leading to the formation of PXR/RXR heterodimer in the nucleus, binding of heterodimer to DNA, mRNA transcription, and translation of mRNA to CYP3A4 proteins (described as a first-order process).
In addition, the model incorporates mRNA background production and degradation, CYP3A4 degradation, and, importantly, a feedback CYP3A4-mediated degradation of rifampicin. First-order kinetics are used to model the diffusion between compartments, heterodimer dissociation, mRNA transcription and translation, degradation of mRNA and CYP3A4, and the formation of the heterodimer, while the degradation of rifampicin by CYP3A4 is modelled by second-order kinetics. Incorporation of zero-order kinetics of mRNA background production and first-order kinetics of CYP3A4 mRNA degradation enables us to model steady-state CYP3A4 mRNA levels. This is very convenient, since there is a constitutive CYP3A4 mRNA as well as CYP3A4 protein expression in normal hepatocytes. In addition, this allows us to use fold-induction data to describe upregulation (induction) of CYP3A4 mRNA or protein from in vitro cellular experiments. These data can be easily obtained using standard qRT-PCR or Western blotting methods. In contrast, absolute quantification of mRNA transcripts is always difficult and needs sophisticated instrumentation or calibration.
The corresponding system of ordinary differential equations for the model is depicted above. Two forms of the model have been proposed in the work. The first one describes the situation in in vitro primary human hepatocytes in culture. It is composed of two compartments: the culture medium outside the plated cell culture (where rifampicin is present) and the interior of the plated cells. The model does not consider the compartmentalization of a cell into a nucleus and cytoplasm, which is a simplification that can lead to error outcomes since mRNA translocation from the nucleus into cytoplasm is an important process in gene regulation.
In the second model, the authors combined the cellular model with a model representing the human body in vivo. The model contains one compartment for blood (with a zero-order intravenous application of rifampicin and taking into account urinary excretion) and one for the liver. This model is parameterized with physiologic parameters, such as blood flow rate, liver and distribution volume, and the liver/blood partition coefficient.
The values necessary to run the models were taken from previously published papers. However, a small number of parameters (4 for the in vitro and 6 for the in vivo model) had to be estimated through least-squares curve fitting. In addition, a sensitivity analysis was performed and computed to assess the influence of the accuracies of the estimated parameters on the output concentrations. All computations were done using MATLAB ® (Math Works ® , Natick, MA, USA) [97].
The models were validated through comparison with published data for rifampicin blood concentrations and CYP3A4 mRNA fold induction data measured in primary human hepatocytes. Whereas the model predictions for rifampicin concentrations in human plasma are acceptable, CYP3A4 mRNA expression profiles at early points modelled in the cellular model are discrepant in comparison with experimental data. Therefore, the authors modified the model and incorporated, through inserting an additional transit compartment, an artificial delay in CYP3A4 mRNA synthesis after PXR activation. In the modified setting, an improved fit of predicted CYP3A4 mRNA expression data with experimental data was observed. The authors suggest that the delay is caused by the fact that the model does not address cytoplasm-to-nucleus mRNA translocation as well as RNA processing, which could be carried out by defining an additional compartment for the nucleus. Moreover, when the dose-dependent induction of CYP3A4 mRNA was modelled, the predictions provided a fine fit with the collected data, but they did not approximate the data accurately at larger doses because a maximum of response was observed due to saturation.
This model thus introduced quantitative modelling of PXR-mediated regulation of the CYP3A4 target gene in a cellular model and enabled rough modelling with different parameters. In addition, the authors tried to scale the cellular model for a liver scenario employing hepatocellularity factor scaling. This attempt, however, was not further applied.

Model by Yamashita et al.
In a paper by Yamashita et al. [88], another mechanism-based model for PXR-mediated human CYP3A4 gene regulation was developed. The paper contains in fact a collection of models for predicting DDIs caused by PXR activation.
First of all, the authors modelled the intracellular induction dynamics of CYP3A4 expression in primary human hepatocytes by three ODEs for the amount of activated PXR, the CYP3A4 mRNA level, and the CYP3A4 enzyme level. As in the previous model, it incorporates mRNA background production and degradation as well as CYP3A4 enzyme-mediated degradation of rifampicin as a negative feedback loop regulation. First-order kinetics are used to model PXR's degradation, to model PXR-induced mRNA transcription and the degradation of mRNA, and to model the mRNA-activated production of CYP3A4 protein and its degradation. The action of rifampicin is not modelled based on diffusion from an outside compartment but using its affinity parameter (EC 50 concentration) with respect to PXR, i.e., using a Hill equation. Zero-order kinetics are used for mRNA background production.
Employing the mathematical modelling and a large set of experimental data from the literature, the authors predicted degradation rate constants for CYP3A4 mRNA and proteins in primary human hepatocytes and showed a mechanistic application of the model in the examination of molecular aspects of PXR-mediated CYP3A4 regulation after a single-dose treatment.
Two extrapolations of this basic model are further presented: An indirect effect dynamic model, which assumes that hepatic CYP3A4 protein levels depend dynamically on the actual liver unbound fraction concentration of rifampicin (yielding just one ODE), and a static model, where hepatic CYP3A4 levels depend on the average steady state concentrations of rifampicin in plasma.
Moreover, a PBPK model for the clinical pharmacokinetics of rifampicin concentration in the blood and the liver is proposed in the paper. Employing the model, the authors simulate the fact that during repeated oral doses, rifampicin clearance is an accelerating nonlinear saturable process due to feedback PXR-mediated induction of rifampicin metabolism. This leads to mass-balance equations not only for rifampicin blood and liver concentrations, but for the inducible maximum elimination rate in the liver as well. The equation for rifampicin liver concentration contains, in addition to the terms in the work by Luke et al. [90], a Michaelis-Menten term and a repeated-dose term.
Computations in the report were performed in the multi-hierarchical physiology simulation platforms PhysioDesigner (for the PBPK models) and CellDesigner (for the intracellular simulation) [98]. They included the estimation of unknown model parameters. To cope with the influence of inter-donor variability of drug metabolism in hepatocytes, least-squares based analyses were carried out using the NONMEM software (Icon plc., Dublin, Ireland) [99].
This model demonstrates that rifampicin-induced drug-drug interactions can be extrapolated with acceptable accuracy using the developed models. Furthermore, the static model can be a tool to simulate DDI once the blood concentration of the inducer reaches steady-state after the administration of repeated oral doses.

Model by Li et al.
A PK/PD model for rat CYP3A orthologues expression in the manuscript by Li et al. [96] describes PXR-mediated induction of CYP3A1/2 enzymes by the rat PXR ligand dexamethasone. It investigates the effect of PXR by monitoring in parallel both CYP3A1 and CYP3A2 mRNAs and proteins expression and by measuring the final total 6β-testosterone hydroxylation formation enzyme activity.
To model dexamethasone concentrations after intraperitonial administration (100 mg/kg) to rats, a two-compartment mammillary PBPK model with two ODEs for the blood and the liver was used. Similarly to the manuscript by Luke et al. (2010), it assumes zero-order drug absorption and comprises systemic (urinal) excretion, but it considers drug concentrations only. The link with PXR activity is indirectly modelled with a Hill equation giving, in dependency of liver dexamethasone levels, the relative concentration of DNA binding sites bound to dexamethasone-PXR complexes. As experimental data in rats suggest that CYP3A1 and CYP3A2 mRNA levels versus time are significantly delayed after dexamethasone application to rats, two parallel chains of transit compartments (filled with abstract substances called stimulation) were used. For CYP3A1, one transit compartment suffices to capture the delay, but CYP3A2 required eight transit compartments. The optimal number of transit compartments was found through stepwise addition or deletion of one transit compartment and was selected from the inflection point on an objective function value (presumably representing a measure of model quality, but this is not specified in the paper) versus compartment number curve.
Next, the stimulation concentration for the final transit compartment predicts mRNA levels using first-order kinetics. CYP3A1/2 mRNA are further assumed to be produced in the background by a zero-order rate constant and mRNA degradation to take place via a first-order process. CYP3A1/2 enzyme activity is then modelled with first-order constants for CYP3A1/2 degradation and synthesis via mRNA translation, where an amplification factor indicates that one copy of the mRNA can be translated into multiple copies of the protein. Finally, the computed protein concentrations of both CYP3A1 and CYP3A2 were combined to predict the resulting overall 6β-testosterone hydroxylation formation.
After fitting pharmacokinetic data, PK parameters were fixed and PD parameters were estimated using NONMEM software version 7.1.2. Applying these models to experimental data, the authors achieved a satisfactory fit of CYP3A1 and CYP3A2 mRNA profiles with parallel estimation of CYP3A1 mRNA and CYP3A2 mRNA delay. Similarly, both CYP3A1 and CYP3A2 protein dynamics with an appropriate lag time between mRNA synthesis and protein synthesis (about 12 h), as well as the CYP3A activity profile, were predicted with sufficient accuracy.
This model very well demonstrates the feasibility of in vivo PXR target gene induction modelling using PK/PD models and might be the most appropriate model that predicts the data rather accurately. This may be due to the fact that the authors used absolute concentration data (but not relative fold-upregulation data) for both CYP3A1/2 mRNA and CYP3A1/2 proteins in the work. The addition of a chain of transit compartments was essential to model the delay between mRNA transcription levels and drug plasma concentrations. Unfortunately, the report is based on animal data and parameters that cannot be easily extrapolated to humans.

Model by Bailey et al.
Recently, Bailey et al. [91] developed a model to describe feedback and feedforward loops of the PXR-mediated transcriptional control of the rat CYP3A1 enzyme involved in steroid hormones deactivation. In the paper, the authors performed experiments with primary rat hepatocytes treated with different concentrations of rodent PXR ligands. They observed expression of whole-genome transcriptome using microarrays and expression of individual target genes using RT-PCR and Western blotting.
The authors formulated two models. The first one describes pregnenalone-16α-carbonitrile (PCN) as a xenobiotic PXR agonist, PXR-induced expression of its target gene CYP3A1, and downregulation of PXR transcript levels. The second model describes the effects of lithocholic acid as a common endogenous agonist for PXR, the farnesoid X-receptor, and the vitamin D receptor (VDR) on regulation of their target genes CYP3A1, Fibrinogen B, and CYP24, respectively. The effect on testosterone, cortisol, and progesterone metabolism was also included in the models. These in silico models were generated using CellDesigner and parameter values were taken from the literature if possible or estimated from experimental data.
The models show that biological response to chemicals occurs by the interaction of NRs via regulatory loops or through competition for ligands and response elements. The levels of PXR, controlled by the feedback loop, are reduced to approximately 50% of that seen with no agonist present. The model demonstrates that this feedback loop is an important factor to prevent hyperexpression of PXR target genes, such as CYP3A, which was confirmed in vitro and may help in explaining the robustness in steroid homeostasis within the cell. Thus, PXR might be regarded not only as a xenosensor but possibly also as an endosensor. The model helps us in understanding how the functions of PXR on xenobiotic and endobiotic metabolism are balanced to achieve both enough sensitivity to eliminate xenobiotics and at the same time maintain the hormonal homeostasis of the organism.

Model by Kolodkin et al.
Kolodkin et al. [89] developed a model to describe the combined regulatory network and interaction of two NR's in the human liver as a response to the stress hormone cortisol, a dual ligand of both PXR and glucocorticoid receptor (GR). As stress exposure causes transient blood cortisol spikes, the authors modelled the effect of cortisol on its two nuclear receptors. Cortisol is considered to be a high-affinity ligand for GR and a low-affinity ligand for PXR. Despite that the model represents only a small part of the whole biological system, excluding the possibility of changes in other proteins outside the model, it gives a realistic and complex view of the regulatory response to glucocorticoids in the human liver, including feedforward and feedback loop regulation.
The modelled GR-induced reactions are upregulation of PXR transcription [100], auto-downregulation of GR [101], and stimulation of the transcription of the classic marker gene tyrosine aminotransferase (TAT) [102], which represents the pharmacodynamic response to GR activation. As for the addressed PXR-induced signals, the formed PXR-ligand complex upregulates a second target gene set, including the marker gene CYP3A4 [103]. The reduction in levels of cortisol by CYP3A4 is modelled in the frame of a feedback loop. Reactions are described using various balance and rate equations. The resulting rather extensive and sophisticated signal regulatory network was constructed by using the CellDesigner software, with parameter values obtained or estimated based on the literature and de novo-generated data.
The model explains, among others, the important role of GR-induced upregulation of PXR. This upregulation allows PXR to bind to a multitude of low-affinity endogenous chemicals, all leading to increased PXR target genes, such as CYP3A4, and as a consequence stimulates the return to homeostasis. Thus, in spite of the low affinity of cortisol for PXR, PXR's promiscuity, in combination with a GR-induced feedforward signal, makes it important for an appropriate response to stress challenge.
The model further describes an essential feature of the network: TAT levels, representing the pharmacodynamic response of GR, are significantly less sensitive to increased frequency of repeated stress episodes than the other modelled species. This demonstrates the ability of the network to adapt and thus limit the magnitude of the physiological response and the risk of disease progression. The model also predicts that there is a negligible difference between the stress response after repeated days with medium stress and the stress response following a weekend of no stress. In contrast, following a short vacation, the model predicts a muted response.

Discussion
Each of the described models has its own merits and restrictions. In this section, we will discuss some options to complement the existing models that address additional mechanistic questions or use different computational methodologies.
The currently published quantitative models of PXR-mediated regulation of its target genes use several approaches. Some models first describe PXR-mediated regulation in cellular models of hepatocytes and then aim to extrapolate the model to an in vivo situation based on indirect correlation of plasma PXR ligand concentrations with compartmental or whole body physiologically based pharmacokinetics models (PKPB) [88][89][90][91]. Some other models use in vivo data obtained in animals [96,104]. Some models are mechanistically based and consider the entry of a ligand into cells or cytoplasm-nuclear translocation of the liganded PXR [90], but others do not consider morphological factors of cells. Importantly, most of the models consider feedback (or feedforward) loops in PXR-mediated regulation of CYP3A genes, but only some of them can implement both of the critical feedback loops, i.e., CYP3A4-mediated degradation of a PXR ligand and PXR-mediated downregulation of PXR gene expression. Interestingly, some models use data obtained after the treatment of hepatocytes in one time interval with gradually increased concentrations of a PXR ligand [89,91].
All models studied the expression of a CYP3A orthologue as the most important target gene of PXR. Of note, Li et al. (2012) found that CYP3A1 and CYP3A2 genes and enzymes are differently regulated by PXR activation in rats, which clearly indicates that PXR target genes can have different time profiles of expression [96].

Feedback Loops
The papers of Kolodkin et al. and Bailey et al. included feedback and feedforward loops into the models. Importantly, PXR autoregulation resulting in PXR-mediated transrepression (downregulation), which was first proposed by Dr. Plant's group [105], has been incorporated into their models.
This autoregulation mechanism has already been applied and demonstrated in the mathematical modelling of quantitative glucocorticoid receptor (GR)-mediated regulation. GR is another member of the nuclear receptor superfamily (NR3C1), which has very similar intracellular transcriptional machinery.
Glucocorticoids are an essential class of drugs indicated for the treatment of autoimmune diseases, such as rheumatoid arthritis. Recently, a fifth-generation model for the corticosteroid dynamic of GR activation has been reported. In the model, which is accompanied by comprehensive animal data, the authors incorporated cellular mechanistic events (including cytoplasm-nucleus recycling of GR) as well as feedback downregulation of GR by glucocorticoids [81]. According to the study, after a single bolus dose of glucocorticoids, the free cytosolic receptor levels fell immediately to zero and returned to baseline by 72 h. When two doses were given at 0 and 24 h, the second dose had a reduced effect due to the downregulation of the GR. Thus, less free cytosolic receptors were available for the next dose. This shows the importance of understanding the downregulation of NRs, including PXR, in multi-dose pharmacotherapy.
According to our knowledge, temporal data concerning the downregulation of PXR was not considered or used in any model. Baily et al. first proposed that PXR negative feedback downregulation of PXR expression attenuates the expression of the CYP3A gene as was shown using in silico modelling and a data set with only 48 h time-endpoint measurements.
We remark that PXR levels exhibit a large inter-donor variability caused by, among other things, differing PXR gene (NR1I2) genotypes, mutations, or variants. We demonstrate the importance of this issue in the validation of the model of Luke et al. with our own experimental data. In the experiment, primary human hepatocytes were treated with 10 µM of rifampicin. The expression level profiles of mRNA for the CYP3A4 enzyme were analyzed using the qRT-PCR method at 0, 6, 12, 24, 48, and 96 h from the beginning of the treatment. The measured average fold levels are displayed as circles in Figure 3 and the blue curves give the levels predicted by the model of Luke et al. [90]. In the left part of the figure (Figure 3A.), we used both the literature parameter values and the estimated parameter values reported by Luke et al. Clearly, our experimental data do not fit well with the predicted estimation. However, one of the estimated parameters is the total amount of PXR in the cell. Doubling the value of this estimated parameter resulted in satisfactory fitting as can be seen in the right part of the figure (Figure 3B.). This example clearly indicates that the mechanistic molecular-based model can be parametrized with individual parameters based on PXR gene (NR1I2) genotype or particular PXR expression. had a reduced effect due to the downregulation of the GR. Thus, less free cytosolic receptors were available for the next dose. This shows the importance of understanding the downregulation of NRs, including PXR, in multi-dose pharmacotherapy. According to our knowledge, temporal data concerning the downregulation of PXR was not considered or used in any model. Baily et al. first proposed that PXR negative feedback downregulation of PXR expression attenuates the expression of the CYP3A gene as was shown using in silico modelling and a data set with only 48 h time-endpoint measurements.
We remark that PXR levels exhibit a large inter-donor variability caused by, among other things, differing PXR gene (NR1I2) genotypes, mutations, or variants. We demonstrate the importance of this issue in the validation of the model of Luke et al. with our own experimental data. In the experiment, primary human hepatocytes were treated with 10 µM of rifampicin. The expression level profiles of mRNA for the CYP3A4 enzyme were analyzed using the qRT-PCR method at 0, 6, 12, 24, 48, and 96 h from the beginning of the treatment. The measured average fold levels are displayed as circles in Figure 3   In the paper by Yamashita et al. [88], the authors used PKPB modelling to analyze the autoinduction metabolic deactivation of rifampicin following repeated oral dosing of rifampicin (clinical data were published in the literature). Thus, the authors show that the feedback inhibition of the PXR ligand rifampicin by CYP3A4-mediated degradation is an important phenomenon that has to be considered in clinics.
Unfortunately, there is little data on the simulation of PXR-mediated feedback loops after repeated dosage during longer time intervals (i.e., at least one week, as has been done for the GR in  In the paper by Yamashita et al. [88], the authors used PKPB modelling to analyze the auto-induction metabolic deactivation of rifampicin following repeated oral dosing of rifampicin (clinical data were published in the literature). Thus, the authors show that the feedback inhibition of the PXR ligand rifampicin by CYP3A4-mediated degradation is an important phenomenon that has to be considered in clinics.
Unfortunately, there is little data on the simulation of PXR-mediated feedback loops after repeated dosage during longer time intervals (i.e., at least one week, as has been done for the GR in Kolodkin et al.). With our Matlab implementation of the in vitro model by Luke et al., we simulated fold-induction of the CYP3A4 enzyme with varying levels of PXR ligand metabolism by CYP3A4 for a period of about two weeks. We perturbed with one order of magnitude the estimated value of 2.47 × 10 −5 in Luke et al. for kmet, a second-order metabolic constant indicating the feedback metabolism of rifampicin in cells. The amount of feedback metabolism seems to strongly influence rifampicin levels only (Figure 4). Predicted CYP3A4 fold-induction levels, especially for small values of kmet, seem very high and will probably not occur in in vivo situations; they might be caused by inaccurate estimation of the precise kmet value or of other parameters of the model. In a similar way, the influence could be predicted of the diffusion speed of the ligand into hepatocytes, of other aspects, or of interactions between them.

Extended Computational Methodologies
A major restriction in mathematical PK/PD modelling is the limited availability of reliable parameters for the model and complications with obtaining data. Some of these complications are connected with technical or methodological difficulties of data measurement in cellular models, in experimental animals, or in humans.
For the estimation of unknown parameters, the described models use curve fitting based on a comparison of observed with simulated concentration time profiles. Mathematically, the fitting consists of the minimization of a sum of squares cost function. Improved accuracy of the estimated parameters can be obtained with a Bayesian approach using Monte Carlo-type simulations. The basic idea is to generate a large number of random values for the model parameters with the help of a suitable pseudo-random number generator. Evaluation of the behavior of the model for random values can provide, using tools from statistics, relatively reliable information on general behavior

Extended Computational Methodologies
A major restriction in mathematical PK/PD modelling is the limited availability of reliable parameters for the model and complications with obtaining data. Some of these complications are connected with technical or methodological difficulties of data measurement in cellular models, in experimental animals, or in humans.
For the estimation of unknown parameters, the described models use curve fitting based on a comparison of observed with simulated concentration time profiles. Mathematically, the fitting consists of the minimization of a sum of squares cost function. Improved accuracy of the estimated parameters can be obtained with a Bayesian approach using Monte Carlo-type simulations. The basic idea is to generate a large number of random values for the model parameters with the help of a suitable pseudo-random number generator. Evaluation of the behavior of the model for random values can provide, using tools from statistics, relatively reliable information on general behavior (assuming that a large number of random values was generated). With the high performance of modern computer systems, such a process often has many advantages over laboratory experiments; it can also complement them or help us in the proper dosing of experiments based on pilot mathematical predictions. An example of a model employing such techniques for the simulation of the dynamic range of 4βHC as a biomarker for CYP3A4 induction and inhibition is presented in the report by Leil et al. [85] and in the paper by Jiang et al. [84].
The models we described are almost completely based on one of the following types of equations: First/zero order kinetics, mass-balance equations, the Hill equation, or Michaelis-Menten kinetics. These equations are considered appropriate for the involved reactions and appear to be used in PK/PD models for most other pharmacological applications as well. Nevertheless, the following remarks might be important. First, the justification for using a particular type of equation is sometimes missing. The process of ligand-binding, for instance, can be modelled with mass-action kinetics if this is known to be the only type of reaction taking place or with a Michaelis-Menten reaction, which is different. The Hill equation can be used to describe the cooperativity of binding to a receptor when more information is known (negative, positive cooperativity). These different modelling options lead to rather different mathematical equations. Different equations for the same portion of the signaling pathway may be due to the fact that some models simply use, from the outset, weaker assumptions than others. However, in some cases, the absence of a biological justification for the used modelling equations might suggest that they have been determined based on a trial-and-error basis; i.e., the equations leading to the best model for the given data were used. Although this seems to be the natural thing to do when no alternatives are available, it would be interesting to verify whether the presented models remain as accurate with different experimental data values for the same phenomena. In other words, the question is whether, in addition to the fitting of unknown parameter values to observed data, one actually fits the modeling equations to the observed data as well. Of course, one can find specific parameters for several models which represent the same given data. However, the goal is to understand and model the phenomena behind and make long-term predictions with the correct model. Here, a wrong model could give a different outcome.
Another example of 'model fitting', with respect to the number of compartments, is the definition of transit compartments to model a transport delay. The number of transit compartments is sometimes determined based on curve fitting to experimental data (e.g., in Li et al. [96]). In fact, introducing transit departments can be very convenient to precisely model reactions that are not understood in detail, i.e., when it is only observed that somewhere, something takes place that delays substance transport without knowledge of the type of reaction or of its location. In the latter case, the transit compartment does not have a precise physical place in the modelled process, which may sound a little counterintuitive. In other situations, complicated delaying processes are willingly replaced by a transit department to keep the overall model simple. If measurements are available, some phenomena in between can be neglected, and a model with limited equations can be better analyzed mathematically. For instance, the nuclear processes of the considered signaling pathway can be neglected because the concentrations leaving the nucleus can be measured. In the model of Luke et al. [90], the encountered delay is explained by the fact that the model ignores nuclear-cytoplasmic compartmentalization to incorporate mRNA nucleus-to-cytoplasm translocation. Hargrove [106] reports that indeed there is a delay before the appearance of mRNA in the cytoplasm that is due to the transcription of the gene, processing of the primary transcript, and nucleocytoplasmic transport [90].
Often, the unexplained or unmodelled reactions in transit compartments correspond, in reality, to the accumulation of a substance at some particular location before it is further transported. Protein synthesis occurs in ribosomes, which are located free in the cytosol or bound to the endoplasmic reticulum where synthesized proteins are localized [107]. Thus, CYP3A4 protein and its mRNA can accumulate near the reticulum. To gain more insight into where reactions take place and what their local effect is, an alternative to defining additional compartments is to incorporate spatial distribution in the substance concentration functions: The functions would be dependent not only on time but space-points as well. This requires an altogether different computational methodology, because differential equations for functions with more than one variable (so-called partial differential equations (PDEs)) have to be considered. Even if substances appear to often be homogeneously distributed, it can be beneficial to provide spatial resolution in compartments where experimental observations suggest heterogeneous distributions. For example, the model of Luke et al. [90] could be modified and the cytoplasm compartment modelled with PDEs (i.e., with spatial resolution) while the nucleus compartment is modelled as usual with ODEs. Mathematically, this leads to a mixed system of PDEs coupled with ODEs. The variables X int (t), PR(t), mRN A(t), and CYP3A4(t) in the equations displayed at the end of Section 5.1 will be dependent not only on a time variable t, but on a space variable x (with possibly three components in three dimensions) as well. For instance, the last equation then takes the form: ∂ t CYP3A4(t, x) = D CYP3A4 ∆CYP3A4(t, x) + k cyp mRN A cyp (t, x) − k cyp,deg CYP3A4(t, x), (6) where D CYP3A4 is the diffusion coefficient (or matrix with anisotropic diffusion) for CYP3A4, ∆ denotes the gradient, and mRN A cyp is the mRNA concentration in the cytoplasm.
PK/PD models including spatial resolution appear to have been proposed only rarely in other applications; some examples are reports by Claus et al. [108] and Friedmann et al. [109]. In the report by Thurley et al. [110], it is shown that high concentration gradients evolve despite fast diffusion. Here, the cytokine concentration is immediately distributed homogeneously in space, but the local reactions control the dynamics and lead to an inhomogeneous steady-state.
Because elevated drug concentrations are often toxic, it can be crucial to monitor not only the average drug level all over a compartment, but to detect possible localized maxima as well. Similarly, the exceeding of no-observed-adverse-effect levels in the context of drug-drug interactions should be detectable locally, inside compartments. The degree of heterogeneity of substance concentrations in an intracellular compartment can depend strongly on the speed of the diffusion, which is often inversely correlated with molecule sizes (e.g., larger molecules, such as rifampicin, penetrate the cell membrane slowly and incompletely). Cell shapes may have an influence as well. In the papers by Claus et al. and Friedman et al. [108,109], the authors found, in a different context, larger cytoplasmic concentration gradients in fibroblast cells than in regularly shaped ones. This can provide a delay or change in the nuclear processes. For these cases, the effect of the cell geometry on the nuclear processes is in fact higher than the effect of diffusion and elaborate diffusion measurements can be relaxed. For modelling signaling transduction in cells with spatial extensions, the more complex PDE modelling should be used. In general, the more complicated the signaling pathway dynamics is, in particular when feedback or feedforward loops are involved, the less obvious it is that homogeneity of concentrations can be assumed: Can a rather small heterogeneity in one of the substances cause a more significant variation in spatial distribution of other substances of interest? Can such a heterogeneity, in conjunction with intricate kinetics, bring about a delay? Even if spatial resolution may not seem necessary at first sight, it might reveal unexpected explanations for observed pharmacological phenomena.

Conclusions
Each of the described models has its own merits and restrictions. Currently, we have several more or less sophisticated models of PXR activation describing the cellular or in vivo scenario of target genes regulation. However, we need additional models based on robust experimental data that would comprehensively describe the molecular, organ, and whole-body aspects of PXR-mediated cytochrome P450 gene expression. In addition, all reports on modelling PXR focused only on the CYP3A4 gene and omitted other important CYP enzymes regulated by PXR. We do not have models that consider additional PXR ligands. Finally, PXR-mediated gene regulation has been exclusively modelled in the liver and the intestine has not been considered.

Conflicts of Interest:
The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.