Abstract
Quantitative mappings were established between drug physicochemical properties (PCPs) and parameter values of a physiologically based, mechanistically realistic, in silico liver (ISL). The ISL plugs together autonomous software objects that represent hepatic components at different scales and levels of detail. Microarchitectural features are represented separately from the mechanisms that influence drug metabolism. The same ISL has been validated against liver perfusion data for sucrose and four cationic drugs: antipyrine, atenolol, labetalol, and diltiazem. Parameters sensitive to drug-specific PCPs were tuned so that ISL outflow profiles from a single ISL matched in situ perfused rat liver outflow profiles of all five compounds. Quantitative relationships were then established between the four sets of drug PCPs and the corresponding four sets of PCP-sensitive, ISL parameter values; those relationships were used to predict PCP-sensitive, ISL parameter values for prazosin and propranolol given only their PCPs. Relationships were established using three different methods: 1) a simple linear correlation method, 2) the fuzzy c-means algorithm, and 3) a simple artificial neural network. Each relationship was used separately to predict ISL parameter values for prazosin and propranolol, given their PCPs. Those values were applied in the ISL used earlier to predict the hepatic disposition details for each drug. Although we had only sparse data available, all predicted disposition profiles were judged reasonable (within a factor of 2 of referent profile data). The order of precision, based on a similarity measure, was 3 > 2 > 1.
The in silico liver (ISL) (Hunt et al., 2006; Yan et al., 2007) (Fig. 1) is a first-generation example of synthetic, discrete, physiologically based analogs that are intended for refining, exploring, and testing hypotheses about the details of hepatic drug disposition. Autonomous components represent spatial aspects of hepatic organization and function. Different drugs can be represented and studied simultaneously or separately. Each ISL component can interact uniquely with any drug-representing object that enters its local environment. The consequences of simulated systemic and local interactions can be measured and studied simultaneously, analogous to how wet-lab experiments are conducted.
In Yan et al. (2008), the simulated hepatic disposition of atenolol, antipyrine, labetalol, and diltiazem, along with coadministered sucrose, used only one parameterized ISL structure for all compounds. A subset of components interacted differently with the particular compounds. Monte Carlo ISL variants simulated compound-specific outflow profiles that matched referent profiles. The results supported two hypotheses. 1) The mappings in Fig. 2 between ISL components and corresponding liver components were sufficiently realistic for the stated model use. 2) The simulated drug-ISL component interaction events mapped to corresponding hepatic disposition events. A goal for this project has been to discover and verify ISL counterparts to the relationships between drug physicochemical properties (PCPs) and pharmacokinetic (PK) mechanisms and use them to make predictions.
A vision motivating research on this class of models is identical to one that has motivated development of traditional physiologically based PK models: by “accounting for the causal basis of the observed data,... the possibility exists for efficient use of limited drug-specific data [italics added] to make reasonably accurate predictions as to the pharmacokinetics of specific compounds, both within and between species, as well as under a variety of conditions” (Rowland et al., 2004). PCP-sensitive, physiologically based PK model parameters necessarily conflate features and properties of the biology (aspects of histology and others) with drug PCPs. In doing so, there is a risk that “the causal basis” becomes obscured because of the conflated biological features that were especially influential in causing some property of the data. Interconnections between sinusoids might be such a feature. In the ISL, because we have built a collection of mechanisms from finer-grained components, we have precise control over conflation, yet the causal basis is still there in the drug-component interaction logic (axioms and rules). Our expectation has been that at some level of granularity, the complexity will be sufficiently unraveled so that the logic for a given ISL component-drug interaction (the “causal basis”) will rely heavily on only a small subset of easily specified, drug and biological attributes. In some cases, information about the molecular attributes may be represented adequately by just one or a few PCPs (or, more specifically, knowledge of the drug-environment interaction for which a PCP is a measure). Once we have validated ISL models that exhibit some of these characteristics, we will have moved closer to the above vision of being able to accurately anticipate the hepatic disposition properties of a new compound, given only its molecular formula and structure (and/or validated PCPs).
Are the mechanisms in the current ISL sufficiently realistic to enable such predictions using only limited data? We hypothesized that the differences in the tuned values of the 10 PCP-sensitive ISL parameters were due primarily to differences in PCPs. If that is the case, then quantitative PCP-to-ISL mappings exist and, when discovered, should enable us to predict the disposition of a new compound given its PCPs. Failure to provide better than ballpark predictions could be taken as evidence invalidating that hypothesis and possibly aspects of the ISL mechanisms. The goal of this project has been to suggest quantitative PCP-to-ISL mappings and then use them to predict the disposition properties of prazosin and propranolol given only their PCPs.
Three different methods were used to implement quantitative PCP-to-ISL mappings: linear regression, an artificial neural network, and a well-established fuzzy clustering algorithm. Each method predicted the PCP-sensitive, ISL parameter values for prazosin and propranolol. Those predicted values were then combined with the already-validated, drug-insensitive ISL parameter values. Simulation of the resulting ISLs gave three different, independent versions of expected hepatic disposition details for prazosin and propranolol along with expected liver perfusion outflow profiles. Those profiles were surprisingly good matches to the observed profiles, strengthening the ongoing validation of ISL mechanisms and supporting the existence of the mapping relationships in Fig. 2. Together they represent an important advance in our ability to predict PK properties for drugs given only their structure.
Materials and Methods
To distinguish clearly in silico compounds and processes from corresponding hepatic structures and processes, we use small caps when referring to the former.
Model Structure and Design. ISL structure is illustrated in Fig. 1. Brief summaries of its novel structure and design follow for convenience. Consult Hunt et al. (2006) and Yan et al. (2008) for additional details. An ISL mimics essential features of the in situ perfused rat liver systems used in Hung et al. (2001) to study the hepatic disposition of the cationic drugs. We assumed that rat hepatic secondary unit function is similar throughout a normal liver and that, on average, the function of different lobules would be indistinguishable.
We therefore represented the entire liver as a collection of lobules. Blood flow is represented as entering lobules via portal vein tracts and vascular septa and draining into branches of a common central vein (CV). The ISL represents each lobule as a network of cell-lined sinusoids connected by solute flows. The network is a directed graph; each node (or vertex) represents a cell-lined sinusoid and a link (or edge) between two nodes represents directional flow. For each lobule simulation, the graph structure is determined stochastically within parameter-specified constraints. Objects representing sinusoidal spaces and functions, along with related spatial features, are called “sinusoidal segments” (SSs). One SS is placed at each graph node. The parameters that control lobule form and function for the experiments in this report are listed in Yan et al. (2008) along with all other ISL features. For convenience, they are also listed in Supplemental Table S1. Only 10 parameters were modified during the course of this study. Except as noted below, all other parameter values were the same as those specified in Supplemental Table S1. The meaning of each parameter is in its name. For example, parameter C2BJumpProb is the probability that a compound, when given the option, will move from Grid C to Grid B; parameter BindersPerCellMax specifies the maximum number of binders that can be contained in any one cell.
A SS is a tube-like structure with a blood “core” surrounded by three identical size spaces in the form of two-dimensional grids. Grid A represents the sinusoid rim. Grid B, which represents endothelial cells and fenestra, is wrapped around Grid A. Grid C, which represents all other spaces, including the space of Disse and hepatocytes, is wrapped around Grid B. Different sinusoids can experience different flows and have different surface-to-volume ratios. Such heterogeneity is represented by specifying two SS classes: S1 and S2. Compared with a S2, a S1 has a shorter internal path length and a smaller surface-to-volume ratio.
We represented the lobule interior as being subdivided into three concentric zones to distinguish the quantitative difference in intralobular structural characteristics and, when needed, enzymatic and transporter gradients. The number of directed graph nodes per zone is always in the order of zone 1 > zone 2 > zone 3. Each lobule zone contains at least one node. Because interconnections between sinusoids are frequent in the periportal region and are rare near the CV, there are more interconnections between zone 1 nodes than between zone 2 nodes. There are no interconnections between zone 3 nodes. All of these features are controlled by parameters (Supplemental Table S1), many of which are stochastic to simulate uncertainty and biological variability.
Simulation of Drug Disposition. ISL experiments follow the same protocol used for in situ perfused rat liver experiments. Specifically, in Hung et al. (2001), a bolus containing similar amounts of the extracellular space marker sucrose and a cationic compound was injected into the portal vein through a catheter (because of this, drug does not enter the liver as a bolus or a short-duration square wave). Consequently, the ISL dosage function is a modified gamma function [given in Supplemental Table S1; see Hunt et al. (2006) for details] rather than an impulse to simulate effects of catheters and large vessels. Perfusate was collected using a fraction collector. The fraction of the administered dose contained per collection fraction was measured. Hepatic compounds of interest in the ISL are represented as compounds (also called drugs): they are objects that move through the lobule. Each compound represents a large number of drug molecules and can interact with each encountered SS feature. A compound's behavior is determined by the PCPs of its referent, along with the consequences of its interaction with the lobule and SS features encountered during its unique trek from portal vein tracts to the CV.
After compounds enter the lobule, they can enter a SS in zone 1 at either the core or grid A. Thereafter, until each is collected at the CV, they have several stochastic options. The parameter CoreFlowRate biases compound movement in the core in the CV direction, simulating blood flow. The parameter SinusoidTurbo biases compound movement in the three spaces in the CV direction, but much less so than CoreFlowRate; it simulates turbulent flow. In the core or grid A, a compound can move within that space, jump from one space to another, or exit the SS. When a compound jumps to grid B, it can move within the space, jump back to grid A, or jump on to grid C. All compounds can move within the extracellular portions of grids B and C. When a compound encounters a cell in grid B or C, it can move into it if allowed to do so by the cell (it is allowed when the parameter isMembraneCrossing is true). When a cell sees sucrose, it also sees that its value of isMembraneCrossing is false; consequently, it will not allow sucrose to enter. All cells in grid C represent hepatocytes. After a compound exits a SS in zone 3, it enters the CV. Its arrival is recorded, simulating being collected by a fraction collector. An animated visualization of these events during a simulation, from start to finish, at both individual ISL and SS levels, for sucrose and anty-pyrine administered together, is available through Supplemental Data of Yan et al. (2008).
The only subcellular functions that were needed to simulate behaviors of the six drugs under study were binding and metabolism. As done in Liu and Hunt (2005, 2006) and Sheikh-Bahaei et al. (2006), transporters can be added when needed. At the start of each simulation, objects functioning as containers and representing cells are placed randomly at some fraction of the available locations in grids B and C. All cellular components that bind or metabolize drugs, such as transporter, enzymes and other cellular materials, are conflated and represented by simple functional objects called binders. Each cell contains a randomly specified (within parameter-controlled limits) number of binders in a well-stirred space. A binder within a grid B cell can only bind and later release a drug. A binder within a grid C hepatocyte, referred to as an enzyme, can bind drug and either release or metabolize it. In the later case, the drug vanishes; for this report, a metabolite replacement is not used.
Because there are 38 ISL parameters (Supplemental Table S1) that specifically influence disposition events, a major simulation task in Yan et al. (2008) was to locate a region of ISL parameter space capable of providing parameter vectors that generate biologically realistic ISL outflow profiles for sucrose and each of the four drugs. To accomplish that, the lobule graph and SS structure was first tuned for sucrose; we held those values constant and then tuned the 10 PCP-sensitive parameters until we matched outflow profiles for each of the four cationic drugs. For this study, the parameters in Table 1 were specified as being sensitive to differences in PCPs. All 28 other ISL parameter values (PVs) are unchanged: they control lobule and simulation features that were considered independent of administered compound and are identical to those specified in Yan et al. (2008). Five of the 10 parameters are used only for compounds that can enter cells. We located a specific ISL parameter domain that could generate valid outflow profiles for all five compounds. The values of the PCP-sensitive parameters provided a base from which we could predict ISL PVs for two new simulated cationic drugs: prazosin and propranolol.
Because of the stochastic nature of ISL simulations, each in silico experiment generates a slightly different outflow profile. In this report, the outflow profiles from 48 simulation trials are summed to represent one experiment.
Mapping PCPs to ISL Parameter Values. Traditional parameter estimation methods (Ljung, 1998) are often infeasible for synthetic, discrete, agent-oriented, biomimetic models of the ISL type, because many of the mathematical assumptions underlying these traditional methods are not satisfied. Therefore, in addition to simply stating and describing the methods used, we also state considerations and provide the rationale.
It is widely held that a mapping exists between the hyperspace of drug PCPs and the hyperspace of measured PK properties (Mager, 2006). Correlation results, as in Hung et al. (2001), frequently demonstrate relationships between a PCP and values of one or more fitted PK parameters, providing evidence of such mappings. Validation of the ISL for five compounds provided evidence for mappings A, B, and C illustrated in Fig. 2: the mapping between ISL outflow profiles and their in situ counterparts (A); ISL events and hepatic disposition events (B); and local ISL mechanisms and their drug-hepatic component counterparts (C). Because a mapping exists between PCPs and the causal hepatic mechanisms, we infer that a logical mapping can also be constructed between drug PCPs and the PCP-sensitive subset of tuned, ISL parameters of the four compounds in Table 1.
How can we establish specific, quantitative mappings so that ISL PVs for a new drug can be approximated directly on the basis of its PCPs? How can we accomplish the task when we have “only limited drug-specific data”? An example of how this can be done for a synthetic, discrete event, agent-oriented model having design features similar to the ISL was recently presented by Sheikh-Bahaei and Hunt (2006) and theirs is one of the methods that we use.
Methods for Predicting ISL Parameter Values. What is the best method and approach for predicting PCP-sensitive ISL PVs for a researcher-selected set of PCPs? We know from exploratory ISL experiments that the change in an outflow profile caused by changing the value of one parameter can be compensated for by modest adjustments in the values of several other parameters: the relationships among ISL PVs and outflow profiles are nonlinear (Hunt et al., 2006). From accumulated PK knowledge, we also know that mappings between PCPs and traditional PK model PVs can also be nonlinear. When one is attempting to identify and unravel complicated mappings, more data are always preferred. Unfortunately, the reality in the PK domain is that because of the high costs of wet-lab experiments and given the combinatorial size of all the factors involved in any PK mechanism, there will always be a shortage of specific data for any particular subsystem. In addition, most PCPs have their own uncertainty issues. Such considerations make predicting ISL parameters challenging. Before suggesting an answer to the lead question, we also need to take into account the following two issues.
-
A goal for this class of models is to incorporate the parameter prediction method(s) into the simulation system so that the ISL and its components can automatically adjust their drug-component interaction rules “on the fly” using the PCP information carried by each compound. This will allow the ISL components to adjust their interactions with drugs even though the particular compound has not been encountered previously, and so no validated interaction rule is available.
-
The ISL is not a static model. It is expected to evolve iteratively with continued use. We can anticipate that components and interaction events will be invalidated and changed. For example, a specific mechanistic detail may be invalidated when the set of targeted attributes is expanded by inclusion of a new attribute (based on wet-lab observations from other experiments). To validate against the expanded set of targeted attributes, a current mechanistic ISL feature may need to be replaced by a revised and possibly more detailed version. Consequently, parameter prediction methods need to be easily transferable to new components as the ISL evolves.
Research in quantitative (drug) structure-PK parameter relationships has considered many mapping methods ranging from simple linear regression to partial least-squares regression [for examples, see Hou et al. (2006) and Ng et al. (2004)] and nonlinear regression methods, including artificial neural networks [e.g., see Nestorov et al. (1998)]. Each has advantages and limitations. For results of the more sophisticated methods to exhibit statistical significance, larger data sets are required, often acquired considering the nature and capabilities of the prediction method. We can use any of these methods as long as we are cognizant of the limitations, especially when working with sparse data: overfitting is an important limitation, but it does not preclude cautious use of a method. Another is the sensitivity of the model's behavior (e.g., the outflow profile) to changes in predicted PVs. The impact of such limitations depends on the intended use to which the predicted outflow profile will be put. Precise, verifiable predictions are beyond the scope of ISL class models [see discussions in Hunt et al. (2006) and Yan et al. (2008)]. At best, an ISL can be used to anticipate the behavior of a new compound in the context of considerable uncertainty, given its assumptions. As an example of the latter: the hypothesized mappings B and C in Fig. 2 are acceptably realistic for the intended use and can be extended to include the new compound. It becomes clear that ISL simulations made using predicted PCPs are a rough approximation of what might be expected if all assumptions are valid. Furthermore, this is not a stand-alone prediction: it should be considered only relative to validated ISL behaviors.
With the full weight of the above caveats in mind, a logical strategy would be to use several parameter prediction methods. Simulations made using each of the predicted parameter sets for the new drug can then be compared with an in situ outflow profile of that same drug. We can state the following. Consider the hepatic disposition properties of prazosin and propranolol determined from wet-lab experiments and analogous properties of the two in an in silico experiment, where the in silico and wet-lab experimental designs are effectively the same. The relative differences and similarities in simulated disposition properties between the two compounds would approximate those of the wet-lab experiment.
Our expectation is that the ISL will be revised and reused to represent the hepatic disposition properties of an increasing number of drugs under a variety of experimental conditions. With improvement, the two hypothesized mappings, B and C in Fig. 2, will become more realistic because validated, finer granularity versions of the simulated mechanisms will be available. As the simulated mechanisms become less abstract, we can expect that a shrinking subset of PCPs (or molecular descriptors calculated from structure information) will become increasingly important in determining the outcome of the referent event. In that case, the differences between predicted, PCP-sensitive PVs, arrived at using different methods, may decrease. If so, the scientific usefulness of the simulations can be expected to increase.
Many molecular descriptors and PCPs (measured and calculated) are available [for examples, see Ng et al. (2004)]. However, motivated by the desire to keep new methods as simple as possible, while still achieving the objective, we elected to focus initially on just the nine in Table 2, the same set used in Yan et al. (2008). We elected not to include pKa values for this study set for the sake of simplicity: we wanted to postpone dealing with the potentially complicated issues of multiple pKa values and the variety of drug-specific pKa values that one can encounter in the literature. In the following, we describe how we used three fundamentally different parameter prediction methods.
Method 1: Correlation and Linear Regression. The simplest method for specifying a quantitative mapping is to use correlation and linear regression. We correlated each of four sets of PCP values in Table 2 against each of four sets of PCP-sensitive, ISL PVs in Table 1. Next, we selected the one correlated pair having the largest correlation coefficient (CC). We then used the corresponding least-squares regression line to predict a value of that ISL parameter for both prazosin and propranolol, given the value of the selected PCP.
Method 2: Fuzzy c-Means Clustering. To handle uncertainty that arises during mappings of this type, we used the fuzzy c-means algorithm (FCMA) introduced by Bezdek et al. (1984) and Pal and Bezdek (1995). It can be used to transform descriptive mappings into quantitative relationships. It starts with influence relationships between subsets of PCPs and subsets of PCP-sensitive, ISL PVs from validated ISLs.
Sheihk-Bahaei and Hunt (2006) used the FCMA and only limited drug-specific data to predict parameter values for their analog of sandwich-cultured hepatocytes. They made the estimation method available to the model to enable it to reasonably anticipate (predict) the biliary transport and excretion properties of a new compound, enkephalin, based on the acceptable, already tuned PVs of salicylate, taurocholate, and methotrexate. By assuming that the similarity among PCPs extends to model PVs, acceptable values of the latter for enkephalin (with PCPs that are very different from the three other drugs) were predicted that had good similarity to measured values. They established a single mapping between 13 selected PCPs and 14 PCP-sensitive model PVs (the ISL has just 10). However, correlations between specific PCPs and specific PK parameter values (Mager, 2006) suggest that a small set of PCPs will have a dominant influence on specific drug-component interactions during disposition, an observation confirmed recently by Ng et al. (2004). It thus seemed reasonable to expect that different subsets of PCPs will map to the hepatic counterpart of different ISL events, such as moving between grids B and C and movement into cells. Furthermore, the artificial neural network method discussed below provides a mapping between all 9 PCPs and all 10 PCP-sensitive, ISL PVs. Therefore, for this second prediction method, we elected to seek usable relationships between subsets of PCPs and the ISL parameters controlling specific ISL events.
Consider the six compounds in Table 2. Each drug's nine PCP values, combined into a vector, can be viewed as a single point in a nine-dimensional hyperspace. The six vectors or points (also referred to by the name of the compound that the PCPs describe) can be classified into one to six groups based in part on their closeness in PCP hyperspace. Any compound can be a member of more than one group: the groups or clusters overlap. The researcher can specify the number of groups or use the FCMA to search for the optimal number. The FCMA assumes that each point has some degree of membership in each group. It then assigns a membership degree value to every vector in each cluster. The FCMA inputs are 1) the data set of n cases to be clustered (for our use, each case is a single drug, a vector of nine PCP values), 2) a parameter m known as the fuzzy exponent, and 3) the number of user specified clusters c. The first two of these inputs are discussed further in the following sections. As recommended by Pal and Bezdek (1995), we always set m = 2. The output of a FCMA, illustrated in Fig. 3, is a c-by-n matrix, containing the membership values of the fuzzy clusters.
ISL Parameter Categories. There are three categories of ISL parameters. Those that control lobule graph and sinusoid structures are the same for all drugs. Those that control the PCP-sensitive parameters comprise the third category. We need only predict the values of PCP-sensitive parameters. We divided the 10 PCP-sensitive ISL parameters in Table 1 into four subgroups: those that reflect the influence of a compound's 1) effective size, 2) movement between spaces, 3) binding to lobular components, and 4) metabolism.
We took the nine PCPs in Table 2 and divided them into five groups, as follows, according to their perceived similarity in influence. Group I includes molecular weight and volume. Group II is the partition coefficient (logPapp). Group III is the unbound fraction (fuB). Group IV is topological parameters: number of rotatable bonds and topological polar surface area. These two PCPs are known to be good descriptors for characterizing drug absorption, including intestinal absorption, bioavailability, Caco-2 permeability, and blood-brain barrier penetration (Ekins et al., 2000; van de Waterbeemd and Gifford, 2003). Group V is hydrogen bond-related parameters: hydrogen bond donor count, hydrogen bond acceptor count, and tautomer count. These three PCPs are known to be important to a compound's chemical interactions with other solutes and migration capabilities in aqueous materials (Ekins et al., 2000; Hansch et al., 2004). Because the relationships between PCPs and ISL parameters are believed to be nonlinear, we conservatively assumed that more than one set of PCPs could influence more than one set of ISL parameters. We then made the following mapping assumptions (Fig. 4) and took them as specifications for moving forward:
-
logPapp is the primary PCP reflecting a drug's partitioning ability and so is expected to map to a compound's ability to move between spaces. We also assumed that parameters in PCP subgroup IV, which reflect drug flexibility and topology, also contribute to a compound's ability to move between spaces: in the ISL, that ability is controlled by ISL parameter subgroup B. We specify that a quantitative mapping can be established between PCP subgroups II and IV and ISL subgroup B. The mapping is illustrated in Fig. 4.
-
In addition to fuB, parameters in PCP subgroup V are expected to contribute to drug binding to (or sequestration in) various cellular components. We therefore specify that a quantitative mapping can be established between PCP subgroups III and V and ISL subgroup C (Fig. 4).
-
drug metabolism (ISL parameter group D) simulates a complicated process and so it may map to a broader range of PCPs. We included the following: logPapp (subgroup II) along with the PCPs in subgroups IV and V.
-
Both molecular weight and volume are known to influence a compound's effective size in biological systems (Hansch et al., 2004). We specify that a quantitative mapping can be established between them and ISL2WetLabScaling.
Fuzzy Cluster Prediction Strategy. In this section, we explain how we used FCMA results to construct quantitative mappings between the PCP and ISL parameter spaces and how we used those mappings to predict the PCP-sensitive, ISL PVs for prazosin and propranolol. For clarity, we first describe the prediction method using an example.
The properties in PCP subgroups II and IV are posited to map to the PCP-sensitive subgroup B parameters (Fig. 4). We illustrate a method to calculate a value of B2CJumpProb for prazosin. The subgroup II and IV values are listed in Table 3. We apply the FCMA to these five three-value parameter vectors, specifying three fuzzy groups: we get the degree of membership in each cluster listed in Table 3. We see that prazosin is primarily a member of a group that contains at least one of the four already-studied drugs. In this case, it is within a group in which labetalol and diltiazem also primary members. Atenolol and antipyrine are primarily in different groups. We assume that the position of prazosin in ISL parameter space relative to the other four drugs is approximately the same as its position in PCP space relative to the other four drugs. As a first approximation, we specify that it is the case. Given that assumption, we can calculate a value of B2CJumpProb using prazosin's PCP degree membership values in each of the three clusters along with the average values of B2CJumpProb for all the members in each of the three clusters. labetalol and diltiazem are members of cluster 1. The average B2CJumpProb value is 0.5; the value for the only primary member of cluster 2 (antipyrine) is 0.35; the value for the only primary member of cluster 3 is 0.3. Given those data, for prazosin we have B2CJumpProb = (0.506)(0.5) + (0.179)(0.35) + (0.316)(0.3) = 0.41. Additional, detailed calculations are provided in the Supplemental Data.
The same approach can be used to calculate the remaining three members of ISL subgroup B. The method can also be extended to calculate expected values for the other ISL parameters. The generalized statement of the method follows.
Parameter Estimation via Analysis of Fuzzy Clusters. In this section, we present an algorithm that uses fuzzy clustering to estimate simulation parameters for a new situation given the tuned parameters of several, previously validated situations. For additional detail, see Sheikh-Bahaei and Hunt (2007). In this work, the new situation is the ISL encountering two new types of objects, one representing prazosin and the other representing propranolol. The algorithm starts with the highest possible number of clusters and iteratively decreases the number of clusters until the new compound is placed in a cluster with at least one of the already-tuned compounds. The compounds in that cluster represent those “most similar” to the new compound.
In general, for a set S, containing n data points S = {c1, c2,... cn} (one data point corresponds to one situation—one drug, in our case, with its associated PCP values), the following algorithm estimates the simulation parameters of a new situation, cn+1:
-
Step 1. Let q = n, and Snew = {c1, c2... cn, cn+1}.
-
Step 2. If q = 1, go to step 4. Else, classify Snew into q clusters using the FCMA.
-
Step 3. If cn+1 is not in the same group with at least another member, then decrease q to q – 1. Repeat steps 2 and 3.
Else, let the G value be the number of group-mates of cn+1.
-
Step 4. Call the q groups G1, G2,... Gq, where cn+1 ∈ G1. Let μk be the membership degree of cn+1 in Gk. Estimate the needed ISL simulation parameters of cn+1 as where PGk is the weighted average parameter vector of all the members of group k:
The accuracy or usefulness of the resulting estimates depends, of course, on how many data points similar to cn+1 exist in the data set, i.e., the higher the G value, the better the accuracy. Because we are starting with just four compounds, four data points, our expectations must necessarily be modest.
Method 3: Artificial Neural Network. To establish a traditional nonlinear mapping, we used a simple feed-forward ANN having a single hidden node (Hudson and Cohen, 1999). The neuron transfer function is y(x) = b/(1 + e–ax); we chose a = b = 1. We recognized that the small training set would lead to an overfitted prediction model. Nevertheless, we used it because the method is known to be good at finding patterns and establishing quantitative relationships between data sets based on those patterns. We also used it to provide an indication of what might be expected from an ANN-based predictive mapping, recognizing that such mappings have often done better than one of the many available multivariate linear mapping options, when the relationships between the data sets are known to be complicated, as is the case here. We used a backpropagation training algorithm that utilized normalized values of the 4 sets of 9 PCPs in Table 2 as inputs; values of the 4 sets of 10 ISL PVs in Table 1 were obtained as outputs. We then used the trained ANN to separately predict a full set of 10 PCP-sensitive, ISL PVs for prazosin and propranolol, given their full set of 9 PCPs as inputs.
Values were normalized as follows: divide each point by its vector Norm. For example, if drug (point) X has PCPs a, b, and c, then it was divided by its length N, i.e., N2 = a2 + b2 + c2 and Xnormalized = (a/N, b/N, c/N). The output values were normalized in the same way. Consequently, the predictions too were normalized; to be used, the predictions had to be scaled: multiplied by the norm of the training vector.
Similarity Measure. The similarity measure (SM) used to compare in silico profiles with the lag-time adjusted, referent outflow profiles is detailed and discussed in Hunt et al. (2006) and Yan et al. (2008). Briefly, we assumed that the coefficients of variation of identically timed observations between t = 1 and t = 100 s from repeat experiments (referent) are constant. For each referent outflow profile measure P, we create two curves, Pl = P(1 – d) and Pu = P(1 + d). They form the lower and upper bounds of a band around P. The SM value is the fraction of evenly spaced, outflow profile measures falling within the band. Here, d is the S.D. of the relative differences between each of six repeated in situ experiments and the mean value for a given collection interval, pooled over all collection intervals. In Yan et al. (2008), in which the ISL outflow profile was iteratively tuned to give a good match to the referent profile, we used d = 1 (S.D.), and an ISL outflow profile was deemed similar to the referent if 80% or more of ISL outflow values were within the band (SM ≥0.8). Here, using predicted PVs without any postprediction refinement, we speculated that having a predicted profile fall within a factor of 2 of the referent profile could be considered a good result. Nevertheless, we calculated SM values using two reference bands: d = 1 and 1.5, which corresponds to factors of 0.45 and 0.675.
ISL Execution. Experiments were executed on an eight-node OSCAR cluster (http://oscar.openclustergroup.org/) running RedHat's Fedora 5. The distribution of the runs uses MPICH 1.2.7 (http://www-unix.mcs.anl.gov/mpi/mpich1/). The ISL was compiled using GCC 4.1.1 against the Swarm 2.2.3 Objective-C libraries (http://swarm.org/wiki/Main_Page). The initial pseudo-random number seed was extracted from the machine's clock. Each completed experiment was archived and logged with the date and time using a Makefile target. Examination and processing of data from simulations used a combination of Matlab (7.0.0) and R (2.4.0) (http://www.r-project.org/) scripts, and Microsoft Excel. As done previously (Hunt et al., 2006; Yan et al., 2008), we smoothed results from each simulation experiment using the Smodulus_smoothing function from the Rwave (version 1.22) package of R. For the results presented, a wavelet window of three observations at 0.5-s intervals was used.
Results
Predicted PCP-Sensitive ISL Parameter Values. Ninety correlations were performed: 4 values of each of 10 ISL PVs (from having already tuned the ISL for four drugs) with the 4 values of each of 9 PCPs corresponding to those four drugs. Table 4 lists the CC values. The regression line for the 10 pairs having the largest CCs (bold) was used to predict the ISL parameter value for both prazosin and propranolol listed in Table 5.
The PCP-sensitive ISL PVs for prazosin and propranolol predicted using eqs. 1 and 2 in combination with the mappings specified in Fig. 4 are also listed in Table 5, along with the ISL PVs for the two drugs predicted using the ANN. Because the referent outflow profiles for prazosin and propranolol included coadministered sucrose, eqs. 1 and 2 were also used to calculate new PCP-sensitive ISL PVs for sucrose. Those values are also listed in Table 5. Several examples of using the FCMA method to predict specific ISL parameter values are provided in the Supplemental Data along with the predicted sucrose outflow profiles (Supplemental Fig. S1) when sucrose was coadministered with prazosin and propranolol.
Hepatic Disposition Predictions. We used the predicted ISL PVs in Table 5 along with the previously specified, PCP-independent, ISL parameter values (Supplemental Table S1) to parameterize new ISLs, three for prazosin and three for propranolol. Simulations using each of these ISLs generated the expected outflow profiles presented in Figs. 5, 6, 7. One simulated outflow profile is the result of pooling 48 separate Monte Carlo trials. One trial represents a single lobule and 48 trials represent the whole liver.
As previously discussed (Yan et al., 2008), the structural and microarchitectural details of the ISL used for each simulation trial were nondeterministic. Even though the number of nodes per zone and the number of intra- and interzone edges were specified, the connectivity pattern was determined randomly at the start of each simulation. SS structures for a given set of parameter values are also stochastic: the actual structures of the 72 SSs used for each ISL run were highly constrained yet different. Because of theses differences, combined with the many forms of discretization, there was variability between ISL instantiations and, consequently, the raw outflow profiles were noisy. That variability reflects, and to some extent simulates, interindividual hepatic differences along with the uncertainty in our knowledge about the exact details of the generative mechanisms. Increasing the number of averaged simulation runs makes outflow profiles smoother. The additional wavelet smoothing was described under Materials and Methods.
As described under Similarity Measure, when matching the simulated outflow profile with lag-time adjusted, in situ outflow profiles, we used a band of ±1 S.D. flanking the original wet-lab outflow profile values as a measure of the similarity of ISL simulations to those wet-lab results. For the outflow profiles generated from the predicted parameter values, two reference bands were used: ±1 S.D. and ±1.5 S.D., and SM values were calculated for each. For these data, a band of ±1 S.D. corresponds to a constant 45% coefficient of variation, and a band of ±1.5 S.D. corresponds to a constant 67.5% coefficient of variation. The ±1 S.D. SM values, in the order prazosin and propranolol, shown in Table 5, were as follows: method 1 (simple linear regression), 0.36 and 0.68; method 2 (the FCMA), 0.50 and 0.58; and method 3 (ANN), 0.79 and 0.97. The corresponding ±1.5 S.D SM values were as follows: method 1, 0.55 and 0.82; method 2, 0.77 and 0.72; and method 3, 0.97 and 0.98.
Drug-Specific Prediction Strategies. Although the predicted ISL parameter values from the FCMA method generated good predicted prazosin and propranolol outflow profiles, they were the result of applying simultaneously a common PCP-to-ISL mapping strategy to both drugs. We recognize that prazosin and propranolol have significantly different PCPs. What if we had sought parameter predictions for each drug separately? Would individualized mappings yield improved outflow profile predictions? The mapping between each compound's PCPs and the matched ISL parameter values are not expected to be identical. We explored the consequences of individualizing the mappings. Several dozen mappings were explored. Starting with the strategy in Fig. 4, we made one mapping change, then predicted parameter values, and simulated outflow profiles. We then repeated the process for a second, third, etc., mapping change. For both drugs, we kept the number of PCP-ISL mappings constant at four, for four subgroups of ISL parameters, as in Fig. 4. We removed or added single PCP values or single PCP groups. Most of the mapping options resulted in expected drug outflow profiles having SM values for both prazosin and propranolol that were smaller than those in Table 5 for the FCMA method. In several cases, the SM values showed little change. Only a few showed improved SM values, but they were modest. The individualized mappings for two of the latter cases (Supplemental Fig. S2) along with the corresponding simulated outflow profiles (Supplemental Fig. S3) are provided in the Supplemental Data.
Discussion
Prediction of PK parameter values plays an increasingly important role in decision making concerning drug candidates in development. The methods and issues, which have been thoroughly reviewed, fall into two categories: empirical, data-based approaches (Stouch et al., 2003; Beresford et al., 2004; Yamashita and Hashida, 2004; Li et al., 2007) and structure-based approaches (Ng et al., 2004; Li et al., 2007). Data-based approaches can be subdivided into linear and nonlinear methods, and each approach may involve clustering. Examples of the former include multiple linear regression, partial least-squares analysis, and linear discriminant analysis. Examples of the latter include ANNs, genetic algorithms, support vector machines, and k-nearest neighbor algorithms. Popular clustering techniques include partitional algorithms (including k-means clustering), nearest neighbor clustering, ANNs, evolutionary approaches, search-based approaches (including simulated annealing), and fuzzy clustering. We used clustering, specifically fuzzy clustering, as one of the three methods because it is difficult to classify PCPs and the drugs they characterize into crisp, disjoint clusters. Conventional clustering produces partitions: each pattern fits in one and only one cluster. Fuzzy clustering assigns each pattern within each cluster a membership value. We speculated that this feature may make it an effective tool for clustering PCPs for mapping to ISL PVs, especially when data are sparse.
The relative merits of linear and nonlinear methods when dealing with traditional PK parameter estimation are now reasonably well understood. However, the ISL is representative of a new class of models. It contains more mechanistic detail than traditional PK models, and only a subset of its parameters is PCP-sensitive. Consequently, we had no a priori basis for selecting one of the preceding methods over another. We chose three quite different methods to begin building insight into their relative merits for parameter estimation within the ISL context, given the multiple levels of uncertainty.
We started with a simple version of each method. We treated the predictions arising from their use as semi-independent “votes” for how to parameterize an ISL, given a set of PCPs of prazosin and propranolol, along with ISL PVs that are independent of PCPs from having validated ISLs for sucrose, altenolol, antipyrine, labetalol, and diltiazem. Simulations arising from use of the three different sets of predicted ISL parameter values gave us different hepatic disposition and outflow predictions. Because referent data were sparse, we were not trying to select a best prediction method. We took each outflow prediction in Figs. 5, 6, 7 as a plausible estimate of what one might expect after administration of prazosin and propranolol to a perfused liver. In the case of the ANN (Fig. 7), we recognized the risks associated with overparameterization (Hudson and Cohen, 1999). All three methods gave rise to predicted profiles that were within a factor of 2 of the actual profiles. We found it noteworthy that the simple, single PCP correlation method (Fig. 5) led to ISL outflow profiles that were close to those of FCMA and ANN.
A single point in discretized PCP space represents a unique compound. The space is not continuous. Each point can map to a relatively unique hepatic outflow profile, which is a consequence of many entangled local, intralobular mechanisms. A mapping from PCP space to a disposition outflow profile (an attribute of the drug-liver phenotype) is complex. A small shift in location in PCP space is expected to correspond to distinguishable outflow profile changes. Table 2 shows that for this work, PCP space has nine dimensions. It is relatively simple to calculate dozens of additional physicochemical descriptors, ranging from simple geometric to sophisticated electrostatic and thermodynamic. On the other hand, with larger data sets, one can reduce the number of dimensions using a linear dimensionality reduction method such as principal component analysis. However, the goal of this work is not finding the optimal size of PCP space. Rather, it is to explore the usefulness of the new methods in this context.
The ISL is an assembly of componentized mechanisms: purposefully separated and abstracted aspects of hepatic form, space, and organization interacting with compounds. Each component mechanism has been unraveled from the complex whole of the hepatic-drug phenotype; it has its own unique phenotype, but that phenotype is much simpler than that of the entire lobule mechanism. Hence, a mapping from component behavior space to PCP space is also simpler and fundamentally different from the PCP-outflow profile mapping described above. For each ISL component, a small location change in PCP space may correspond to a negligible or modest change in the properties of that component mechanism, as well as to the value of the parameter controlling its behavior. We see evidence of this in Tables 1 and 5: there are cases in which different compounds have similar ISL PVs. Note also that for both prazosin and propranolol, there are cases for which each of the three methods predicted similar PVs for the same drug. Therefore, considered in retrospect, it is not surprising that even an overly simple estimation method can lead to an acceptable parameterization. Interestingly, the current level of mechanistic granularity is sufficiently fine-grained so that the subtle patterns in PCPs and in the PCP-sensitive, ISL PVs that are needed to achieve larger SM values (Fig. 7; Table 5) are just enough for the simple ANN parameter estimation technique, the more complicated of the three, to find them.
The results suggest that when use of the synthetic method of assembling separated aspects of mechanism, a parameter estimation method, which reasonably quantifies the relative differences between compound-specific behaviors at the level of detail represented by those mechanisms, will provide satisficing (the minimum satisfactory outcome) ballpark estimates. That finding bodes well for using the ISL class of models for predicting PK properties, given only molecular structure information. However, the six drugs in this study are all weak bases. As such, they occupy a common subregion of PCP space. More drug cases, spanning additional regions, will be needed to build confidence that the above apparently favorable parameter prediction situation is a consequence of the synthetic method and its instantiation for the ISL (rather than fortuity).
Because it is a synthetic model, the ISL successfully shrinks the PCP-sensitive phenomenal manifold relative to that of tightly coupled mathematical PK models. Traditional equation-based, PK models reduce the phenomenal manifold to a set of real-valued parameters interacting within and assuming a relational continuum. Although one can reduce the number of parameters being considered with these models, to maintain accurately quantified relationships, the set of reduced parameters must still accurately capture the aggregated relationships of the larger parameter set. This implies that, for tightly coupled PK models, the phenomenal manifold remains just as complex even though there are fewer parameters in the PK equations. The parsimony achieved by such reduction is, in a sense, illusory. The fine-grained relationships that traditional PK models conflate into mathematical parameters are deeply embedded in the coarse-grained relationships of these parameters. In contrast, the ISL can be actually simplified because it is an assemblage of coupled, abstract mechanisms, such as transfer between spaces and entry into cells, wherein the fine-grained relationships can be completely ignored.
Of course, a negative consequence of the ISL's componentized simplification is that the model can sometimes yield dramatically inaccurate and abiotic simulations. Here are three examples: 1) A parameterization makes too much space available to drugs causing a prolonged, abiotic lag-time to occur before a slowly rising wave of drug emerges. 2) When the diversity in travel paths becomes too limited, travel paths can form distinct clusters, rather than being smoothly distributed. That discontinuity in path options can cause abiotic oscillations to occur in outflow profiles. Evidence of oscillations can be seen in the propranolol profile prediction in Fig. 6. 3) Inaccurate and abiotic crowding can occur when too many drugs converge too quickly on a too small space; this can occur when too many edges connect to one SS. However, all models are inaccurate. A central purpose of the synthetic method is to clearly identify and denote both the accurate and inaccurate assumptions and behaviors that characterize a model. Similar phenomena occur with all wet-lab models. Aspects of the in vitro phenotype are expected to map to corresponding attributes in animal species or patients. However, combinations of experimental specifications can cause in vitro behaviors and phenomena for which there is no counterpart in the referent.
As long as we maintain awareness of the above negative consequence, such effects can be minimized. For research that adopts the motivating vision of physiologically based PK, it seems important to use both the synthetic method for chronicling modeling assumptions (heuristic value) and the inductive method to achieve a high-tolerance predictive value. We posit that, by having simplified, yet fine-grained, heuristic models, such as the ISL, which are commensurate with relatively simple parameter estimation methods, the biological mechanisms become more understandable and accessible. Increased accessibility can provide the methodological leverage needed to build synthetic analogs that will enable us to explore, digest, categorize, understand, and use the massive amounts of data being generated by “omic” technologies.
Acknowledgments
We thank our collaborator Michael S. Roberts for providing us with the original liver perfusion data; we also thank members of the BioSystems Group for helpful discussion and commentary.
Footnotes
-
This work was supported in part by grants (to C.A.H.) and fellowships (to S.S.-B. and S.P.) provided by the CDH Research Foundation and abstracted in part from the Ph.D. dissertation presented by L.Y. to the Graduate Division, University of California, Berkeley, CA.
-
Article, publication date, and citation information can be found at http://dmd.aspetjournals.org.
-
doi:10.1124/dmd.107.019067.
-
ABBREVIATIONS: ISL, in silico liver; PCP, physicochemical property; PK, pharmacokinetic; CV, central vein; SS, sinusoidal segment; PV, parameter value; CC, correlation coefficient; FCMA, fuzzy c-means algorithm; ANN, artificial neural network; SM, similarity measure.
-
↵ The online version of this article (available at http://dmd.aspetjournals.org) contains supplemental material.
- Received October 28, 2007.
- Accepted January 24, 2008.
- The American Society for Pharmacology and Experimental Therapeutics