This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.104.045633v1 88/1/715    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Tran, T. T. Articles by Bentz, J. Search for Related Content PubMed PubMed Citation Articles by Tran, T. T. Articles by Bentz, J.

The Elementary Mass Action Rate Constants of P-gp Transport for a Confluent Monolayer of MDCKII-hMDR1 Cells

Thuy Thanh Tran ^* , Aditya Mittal ^*, Tanya Aldinger , Joseph W. Polli , Andrew Ayrton , Harma Ellens  and Joe Bentz ^*

^* Department of Bioscience & Biotechnology, Drexel University, Philadelphia, Pennsylvania 19104; Preclinical Drug Metabolism and Pharmacokinetics, GlaxoSmithKline, King of Prussia, Pennsylvania; Preclinical Drug Metabolism and Pharmacokinetics, GlaxoSmithKline, Research Triangle Park, North Carolina; and Preclinical Drug Metabolism and Pharmacokinetics, GlaxoSmithKline, Welwyn, England

Correspondence: Address reprint requests to Dr. Joe Bentz, Fax: 215-895-1273; E-mail: bentzj{at}drexel.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS KINETIC MODEL OF TRANSPORT... ALGORITHM FOR HIERARCHICAL... RESULTS DISCUSSION APPENDIX: GLOSSARY ACKNOWLEDGEMENTS REFERENCES

 
The human multi-drug resistance membrane transporter, P-glycoprotein, or P-gp, has been extensively studied due to its importance to human health and disease. Thus far, the kinetic analysis of P-gp transport has been limited to steady-state Michaelis-Menten approaches or to compartmental models, neither of which can prove molecular mechanisms. Determination of the elementary kinetic rate constants of transport will be essential to understanding how P-gp works. The experimental system we use is a confluent monolayer of MDCKII-hMDR1 cells that overexpress P-gp. It is a physiologically relevant model system, and transport is measured without biochemical manipulations of P-gp. The Michaelis-Menten mass action reaction is used to model P-gp transport. Without imposing the steady-state assumptions, this reaction depends upon several parameters that must be simultaneously fitted. An exhaustive fitting of transport data to find all possible parameter vectors that best fit the data was accomplished with a reasonable computation time using a hierarchical algorithm. For three P-gp substrates (amprenavir, loperamide, and quinidine), we have successfully fitted the elementary rate constants, i.e., drug association to P-gp from the apical membrane inner monolayer, drug dissociation back into the apical membrane inner monolayer, and drug efflux from P-gp into the apical chamber, as well as the density of efflux active P-gp. All three drugs had overlapping ranges for the efflux active P-gp, which was a benchmark for the validity of the fitting process. One novel finding was that the association to P-gp appears to be rate-limited solely by drug lateral diffusion within the inner monolayer of the plasma membrane for all three drugs. This would be expected if P-gp structure were open to the lipids of the apical membrane inner monolayer, as has been suggested by recent structural studies. The fitted kinetic parameters show how P-gp efflux of a wide range of xenobiotics has been maximized.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS KINETIC MODEL OF TRANSPORT... ALGORITHM FOR HIERARCHICAL... RESULTS DISCUSSION APPENDIX: GLOSSARY ACKNOWLEDGEMENTS REFERENCES

 
Multi-drug resistance transporters are clinically important because of their ability to render cells resistant to many chemotherapeutic agents, such as anticancer drugs and antibiotics (Borst and Elferink, 2002; Gottesman, 2002; Poelarends et al., 2002; Ambudkar et al., 2003). In addition, it has become increasingly clear over the past decade that some mammalian multi-drug resistance transporters also play an important role in the absorption, distribution, and elimination of drugs and xenobiotics, and may be responsible for clinically important drug-drug interactions (Lown et al., 1997; Schinkel, 1998; Goh et al., 2002). Binding and transport of substrates by P-gp (the P-glycoprotein product of the human MDR1 gene) has been intensely studied, providing important insights into the structure and function of the protein (Gottesmann and Pastan, 1993; Senior et al., 1995; Ambudkar et al., 1999, 2003; Loo and Clarke, 1999; Hrycyna, 2001; Borst and Elferink, 2002; Schmitt and Tampe, 2002; Seelig and Gatlik-Landwojtowicz, 2004). The catalytic cycle of P-gp and its conformational changes have been studied using purified P-gp remixed with lipids or reconstituted proteoliposomes (Sharom et al., 1993; al-Shawi et al., 2003; Loo et al., 2003; Urbatsch et al., 2003). There is a dose-dependent response between ATPase activity and transport activity in cells (Ambudkar et al., 1997) and reconstituted P-gp proteoliposomes (Sharom et al., 1993; Omote and al-Shawi, 2002).

Any model of the P-gp catalytic cycle must be tested against a rigorous kinetic analysis of transport activity for P-gp expressed in a physiologically relevant cell system and, eventually, in vivo. To date, kinetics analyses of P-gp transport have been limited to steady-state Michaelis-Menten approaches (Stein, 1997; Ho et al., 1995; Doppenschmitt et al., 1999; Kuh et al., 2000; al-Shawi et al., 2003; or compartmental analyses (Demant et al., 1990; Ashida et al., 1998; Ito et al., 1999), neither of which can prove transport mechanism.

Our experimental system was the confluent monolayer of MDCKII-hMDR1 cells that overexpress P-gp (Evers et al., 2000; Tang et al., 2002a,b; Tran et al., 2004). This system was chosen because it is a physiologically relevant model for P-gp in tissue, and it involves no biochemical manipulations to measure transport. Substrate can be added to either face of the polarized cell monolayer, apical or basolateral, and transport to the other chamber is monitored. A potent P-gp inhibitor was used to independently measure passive permeability through the cell monolayer. The P-gp substrates used were amprenavir, an HIV protease inhibitor; quinidine, a Na⁺ channel blocker; and loperamide, an antidiarrheal drug. These drugs were chosen because they are good P-gp substrates and show different mass balance problems (Tran et al., 2004). Furthermore, amprenavir shows no evidence for transport saturation (Polli et al., 2001; Bentz et al., unpublished data), rendering the standard Michaelis-Menten analysis of P-gp transport useless. So amprenavir provided a challenging test case.

A minimal comprehensive mass action model for P-gp transport through a confluent cell monolayer must include the known kinetic barriers to transport. This means that several parameters must be fitted simultaneously, and each has a wide range of potential values. This generates a very large ensemble of initial parameter guesses to be fitted against data. Thus, a hierarchical algorithm for fitting was required to get fits in a reasonable time. This computationally intensive approach is the only way to analyze the kinetics of a system with three or more important elementary parameters. We have applied this approach successfully to the kinetics of membrane fusion (Bentz, 2000; Mittal and Bentz, 2001; Mittal et al., 2002; Bentz and Mittal., 2003). Nevertheless, the first hypothesis of this study was that the data from the experimental system would be consistent enough and the kinetic analysis robust enough that a small compact set of parameter values yielding best fits could be found at all.

The mass action model for P-gp transport is precisely the same as used by all prior studies, i.e., the Michaelis-Menten reaction. Our analysis is more complex because we fit all of the parameters of this reaction. Past work has just used the Michaelis-Menten (or Eadie-Hofstee or Wolf-Hanes) steady-state equations to fit the data. For some data, e.g., the ATPase activity of P-gp, the steady-state equations appear to work well enough. The ATPase of P-gp is in the water, like a soluble enzyme. For other data, e.g., transport through a confluent monolayer of cells, we have found that these equations yield very inaccurate predictions for the elementary parameters (Bentz et al., unpublished data). This is why we focus only on the elementary parameters here.

There is a tendency to view a more complicated analysis of a well-known problem as a way of curve fitting, with enough parameters to make the fit look better. However, aside from the rate constants of the Michaelis-Menten reaction and the density of efflux active P-gp, every other parameter we introduce was fitted to independent data with the minimum flexibility possible. The rate constants of the Michaelis-Menten reaction and the density of efflux active P-gp had to be fitted simultaneously, just as would be the case for a soluble enzyme when the active enzyme concentration is not known. To impose the steady-state equations onto such data would be a choice based upon assumptions. A rigorous fitting is always better and is now possible.

Our findings suggest new hypotheses and, at the same time, support the simplest picture for a transporter whose job is to remove a wide range of amphipathic xenobiotic compounds from epithelial cells as rapidly as possible. First, it seems clear that the primary selective pressure on P-gp function was to maintain a wide substrate range, which means that the binding constants for nearly all substrates should be weak. This is what we found for all three drugs. In fact, were it not for the substantial partition coefficients of these drugs into lipid bilayers, their membrane concentrations would be too low to show significant efflux. Second, for P-gp to acquire the substrate rapidly, the on-rate to P-gp from the inner apical monolayer should be fast. All three drugs had very large association rate constants, suggesting that lateral diffusion of the drug through the lipid bilayer was the sole rate-limiting step. Third, the P-gp efflux rate constants for quinidine and loperamide were about the same as the currently known maximal rate of ATP hydrolysis by P-gp, whereas that for amprenavir was substantially larger. We also note that the combination of weak binding to P-gp and ATP hydrolysis-limited efflux requires that efflux of a substrate into the apical chamber occurs only after many tens of thousands or even millions of substrate molecules have visited the binding site and returned to the inner apical monolayer.

The fitted value of the density of efflux active P-gp, was substantially smaller than values measured by quantitative Western blots in similar, but not identical, cell lines. We believe the reason is not due to misfolded P-gp or substantial intracellular sequestration, but rather due to the architecture of the microvilli of the apical membranes of these cells. We speculate that efflux into the apical chamber happens only from the P-gp at the tips of the microvilli. Drug exiting from P-gp at the base of a microvillus would undergo a long random walk, involving many encounters with the same or a nearby microvillus and the P-gp therein. The quantitative Western blot gives whole cell P-gp estimates, whereas our kinetic analysis estimates only the amount of P-gp that actually succeeds in effluxing the drug directly into the apical chamber. This hypothesis is being studied currently.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS KINETIC MODEL OF TRANSPORT... ALGORITHM FOR HIERARCHICAL... RESULTS DISCUSSION APPENDIX: GLOSSARY ACKNOWLEDGEMENTS REFERENCES

 
Materials
Amprenavir and GF120918 were from GlaxoSmithKline, USA (Research Triangle Park, North Carolina). Loperamide was from Sigma-Aldrich (St. Louis, MO) and quinidine from Fisher Scientific (Hampton, NH). ³H-Loperamide (10 Ci/mmole) and ³H-amprenavir (21 Ci/mmole) were custom-synthesized by Amersham Pharmacia Biotech (Little Chalfont, England). ³H-Quinidine (20 Ci/mmole) was from ICN Biomedical (Costa Mesa, CA). Dimethylsulfoxide was from Sigma-Aldrich. Dulbecco's modified Eagle's medium (DMEM) was from MediaTech, VWR (Herndon, VA). Transport medium (DMEM with 25 mM HEPES buffer, high glucose, L-glutamine, pyridoxine hydrochloride, without sodium pyruvate, and without phenol red) was from Gibco (Carlsbad, CA). Cholesterol and porcine brain lipids were from Avanti Polar Lipids (Alabaster, AL). Transwell 12-well plates with polycarbonate inserts were obtained from Costar, (Acton, MA). Ultima Gold scintillation cocktail was from PerkinElmer Life Sciences (Boston, MA).

Cell line and culture conditions
The Madin-Darby canine kidney cell line, which overexpresses human MDR1 (MDCKII-hMDR1), was purchased from the Netherlands Cancer Institute (Amsterdam, The Netherlands). Cells were split twice a week and maintained in culture medium (DMEM supplemented with 10% fetal bovine serum, 50 units/ml penicillin and 50 µg/ml streptomycin). Cells were kept at 37°C in 5% CO₂.

Transport assays
Cells were seeded in 12 well plates at a density of 200,000 cells per insert and grown for 4 days in culture medium. Preliminary results using higher and lower plating densities suggested that this density was optimal with respect to obtaining transport stable confluent monolayers at 4 days. After many preliminary experiments, the data presented here for amprenavir were taken from cell passages No. 32 and No. 33 (50–90 µM), No. 34 (100 µM) and No. 36 (1 and 10 µM). For quinidine, No. 36 and No. 37 were used, and for loperamide, No. 39 and No. 40 were used. Cells were given fresh media 1 day after seeding.

Before the experiment, culture medium was removed and cells were preincubated for 30 min with either transport medium alone (see above) or transport medium supplemented with 2 µM GF120918 to inhibit P-gp. Transport of a range of concentrations of amprenavir, loperamide, and quinidine across the confluent monolayer of cells was measured in both directions, i.e., apical to basolateral (A > B) and basolateral to apical (B > A) in the presence and absence of GF120918. For incubations in the presence of GF120918, the inhibitor was added to both chambers; 0.5 µCi/ml of ³H-amprenavir, ³H-quinidine, or ³H-loperamide was added to each respective drug concentration to allow quantification of transport from donor to receiver chambers by liquid scintillation counting. In addition, ¹⁴C-mannitol (0.75 µCi/ml) was added to monitor cell monolayer integrity. At the indicated time points, 25 µL samples were taken from both donor and receiver chambers, mixed with 10 ml of Ultima Gold scintillation cocktail and counted using a Hewlett Packard Liquid Scintillation Counter. The first time point taken was after 6 min, and we used these data at the starting point for fitting. This eliminated the need to correct for some drug binding to the Costar Transwell apparatus and other initial transient effects (Tran et al., 2004). All simulations accounted for the 25 µL aliquots, since loss of volume was not negligible, especially in the apical chamber (Tran et al., 2004).

Cell stability and substrate metabolism
Tran et al. (2004) showed that the stability of the cell monolayer and plasma membrane with respect to passive and active transport was not affected by the prolonged exposure times to amprenavir for at least 6 h. It was also shown that metabolism or decomposition of substrates was insignificant for all drugs on this timescale using Radio HPLC (Agilent Technologies, Palo Alto, CA).

Liposome preparation
Three liposome preparations were used to roughly mimic the lipid compositions of the relevant monolayers of the plasma membrane (Hill and Zeidel, 2000), while not exceeding ternary mixtures:

1. The apical membrane outer monolayer mimic was a (1:1:1) mixture of phosphotidylcholine/sphingomyelin/cholesterol, denoted PC/SM/chol.
   
2. The basolateral membrane outer monolayer mimic was a (2:1) mixture of phosphotidylcholine/cholesterol, denoted PC/chol.
   
3. The plasma membrane monolayer, facing the cytosol, mimic was a (1:1:1) mixture of phosphotidylserine/phosphatidylethanolamine/cholesterol, denoted PS/PE/chol.
   

Lipid mixtures were dried into thin lipid films in round bottom flasks (Buchi Rotavapor R-200, Brinkmann Instruments, Herisau, Switzerland), and freeze-dried overnight on a Flexi-Dry MP freeze-dryer (Kinetics, Union City, CA). Lipid films were suspended in phosphate buffered saline pH 7.4 (Gibco), freeze-thawed in methanol over dry ice for 20 cycles and heated to 35–40°C, and stored at 4°C overnight. Extrusion was done with a Lipex thermal extruder (Northern Lipids, Vancouver, Canada) through two 0.08 µm Nucleopore membranes for at least seven cycles at room temperature. For the PC/SM/chol liposomes, the first three cycles were run at 65°C. Sizes of liposomes were measured by a NICOMP 380-Submicron Particle Sizer (NICOMP, Santa Barbara, CA). The average diameters for all compositions were 100 ± 30 nm. The amount of lipid recovered after extrusion was determined by the absorbance reading of rhodamine-PE, added to each lipid mixture at 0.1 mol %. Lipid recovery as measured by rhodamine-PE was 80–82%.

Equilibrium binding study
Binding of drugs to liposomes was determined by equilibrium dialysis using a 20-cell equilibrium dialyzer (Spectrum, Fort Lauderdale, FL) and the Spectra/Por 4 membrane (Spectrum), which has a 12–14 kD molecular mass cutoff. The partitioning of 10 µM of each drug (amprenavir, loperamide, and quinidine) into 10 mM liposomes was examined. Each half of the Teflon cells received equimolar concentrations of the cold drug. The half cell with liposomes received an additional 0.25 µCi/ml of radiolabeled drug. The cells were allowed to equilibrate in a 37°C water bath, and samples were taken from both half cells at 8 h. Preliminary studies showed that partitioning was independent of drug concentration (for the ranges used in this work), lipid concentration (5–10 mM) and time after 6 h. At the indicated time points, 25 µL samples were removed from the donor and receiver chambers, mixed with 10 ml of Ultima Gold scintillation cocktail and analyzed by a Packard TriCarb 3100-RT liquid scintillation counter.

Partition coefficients were calculated as K = Cm/Cw, where Cw was the mol of drug per liter of aqueous buffer and Cm was the mol of drug per liter of lipid using the average specific volume of 1.6 µL/µmol of total lipid (Chen and Rand, 1997). The equilibrium partition coefficients measured are shown (see Table 2). The molecular basis for differences is under investigation. However, they do show that for association time constants of seconds, the time constants for release from the liposomes is on the order of minutes, roughly validating our assumption of the partition equilibrium over the hourly time course of transport.

	KINETIC MODEL OF TRANSPORT ACROSS A CONFLUENT CELL MONOLAYER

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS KINETIC MODEL OF TRANSPORT... ALGORITHM FOR HIERARCHICAL... RESULTS DISCUSSION APPENDIX: GLOSSARY ACKNOWLEDGEMENTS REFERENCES

 
MDCKII-hMDR1 cells polarize with the basolateral membrane attached to the polycarbonate filters (Butor and Davoust, 1992; Tran et al., 2004). Fig. 1 is a cartoon of a confluent cell monolayer, with P-gp (upward arrows) expressed on the apical surface. The apical and basolateral chambers are kept separate by the tight junctions. Van Meer and Simons (1986) showed that the tight junctions of the MDCK cell monolayer kept the apical membrane outer monolayer lipids separate from the outer basolateral membrane monolayer lipids while there was significant lateral diffusion throughout the inner plasma membrane, i.e., between the apical and basolateral sides. Our results quantitatively confirm their results.

View larger version (30K): [in this window] [in a new window]   FIGURE 1  Model of a confluent cell monolayer, with the apical membrane on top and the basolateral membrane on the bottom, where it binds to the polycarbonate insert. Passive permeability occurs in both directions. P-gp expressed on the apical membrane transports drug from the inner apical membrane monolayer into the apical chamber. There are two kinetic barriers to transport (see Fig. 3 and accompanying text) flip-flop across the basolateral and apical membranes, with P-gp catalyzing the latter step. The concentration of drug in the apical and basolateral chambers, CA and CB, are measured, while the concentration of drug in the inner plasma membrane, CPC, and the cytosol, CC, are predicted as part of the data-fitting process.

View larger version (17K): [in this window] [in a new window]   FIGURE 3  Passive and active transport of amprenavir across the monolayer of confluent MDCK cells is shown by nmol in the receiver chamber over 6 h when the donor side begins with 90 µM amprenavir. The solid symbols show passive permeability across the cell monolayer when the basolateral chamber was the donor side B > A () and when the apical side was the donor A > B (•). The active transport by P-gp was completely inhibited by preincubation with 2 µM GF120918 added to both chambers. The open symbols show the effect of P-gp on transport of amprenavir across the cell monolayer, which occurs in the absence of GF120918. The fact that the activation of transport in the B > A direction () over the passive transport by P-gp, indicated by the up arrow, is greater than the inhibition in the A > B direction (O), indicated by the down arrow, proves that there are only two kinetic barriers to transport, as explained in the text.

Our aim is to construct a kinetic model for the overall transport process, which contains the minimum number of essential steps to adequately characterize the whole process at a molecular level. The novelty of this model lies primarily in having P-gp as an explicit chemical species, rather than embedded within the Vmax of the Michaelis-Menten analysis or represented by a "rate constant" between compartments. This makes the surface density of efflux active P-gp an explicit variable of the model, to be fitted exhaustively along with the elementary rate constants. Since the fitted range for the surface density of efflux active P-gp should be the same for all drugs, we have a concrete benchmark for the validity of the mass action model and the fitting algorithm.

Another novel element of our analysis was the method for extracting the transport due to P-gp from the total transport, i.e., passive permeability plus active transport. We explicitly fit passive permeability coefficients when P-gp is fully inhibited and then use those rate constants in the fitting of the total transport, when P-gp is active. All prior work has simply subtracted the passive transport data curve from the total transport data curve, and defined this difference as that due to P-gp active transport alone. This would be correct only if transport were irreversible, which is certainly not the case for our experimental system. For example, when there is active transport and the basolateral chamber is the donor, then at any time there is a greater concentration of substrate in the apical chamber than with passive permeation alone. Therefore, there is greater backflow back to the basolateral chamber. The subtraction of passive data from total data underestimates the P-gp transport because backflow is ignored.

Because of these novel approaches in both passive and active kinetics, a fairly detailed description is needed here. First the passive transport model will be constructed, followed by the active transport model. The bottom line for the analysis of passive permeability is to accurately remove this transport from the total transport, so that P-gp mediated transport is accurately estimated. To keep the terminology straight, a glossary is included in the Appendix.

The passive permeation kinetic model
In Tran et al. (2004), we derived an exact solution for the passive permeation across a single permeation barrier, and used this equation to fit the passive permeability coefficient for drugs across the confluent MDCKII-hMDR1 cell monolayer in the presence of the P-gp inhibitor. That is, we measured the overall passive permeability coefficient for permeation from the basolateral chamber to the apical chamber, denoted P_BA, and for the opposite direction, denoted P_AB. We included the possibility of substrate loss by a first order process, e.g., binding to receptors, the experimental apparatus, hydrolysis and/or metabolism. Although the cell monolayer clearly has at least two permeation barriers, the apical and basolateral membranes, the characterization as through a single barrier is adequate for most applications, since intracellular concentrations are not measured. Here we must expand the treatment of passive permeability to include explicitly the apical and basolateral membranes, since P-gp picks up drug from the apical membrane inner monolayer (Sharom et al., 2001; al-Shawi et al., 2003).

Our major assumption was that the substrate concentration in an aqueous compartment is at equilibrium with the lipid monolayer facing the compartment, i.e., we can use a partition coefficient to relate the two concentrations. This is a reasonable assumption for this system, where transport occurs on an hourly timescale. For example, Abreu et al. (2003) have shown that the fluorophor Rhodamine Green-tetradecylamide associates with lipid bilayers within a few seconds. The assumption of partition equilibrium reduces the number of equations to be fitted by about a factor of two. This assumption also means that drug recycling is maximal, since drug effluxed into the apical chamber is immediately equilibrated with the outer apical monolayer.

Equilibrium dialysis binding of drugs to intact cells cannot estimate these individual cell membrane monolayer partition coefficients very well because these drugs are all permeable enough to bind throughout the cell within an hour. Therefore, we used 0.1 µm extruded unilamellar liposomes (LUV) whose lipid compositions mimic, in a very simple way, the lipid compositions of the respective membrane monolayers. For example, we use the partition coefficient for the LUV composed of sphingomyelin/phosphatidylcholine/cholesterol (SM/PC/chol, 1:1:1) to relate the concentration of substrate in the apical chamber, C_A, which we can measure, to the concentration in the apical membrane outer monolayer, C_AO, which we can't measure. We use PC/chol (2:1) as a mimic for the outer basolateral monolayer and phosphatidylethanolamine/phosphatidylserine/chol (PE/PS/chol, 1:1:1) as a mimic for the inner plasma monolayer, i.e.,

(1)

K_BO, K_AO, and K_PC denote the respective partition coefficients, while the LUV lipid compositions used are given in parentheses. Obviously, cell membranes have much more complex compositions (Hill and Zeidel, 2000), but we did not want to exceed ternary mixtures or use exotic lipids.

The moles of drug in contact with the basolateral side of the B > C kinetic barrier are V_BC_B + V_BOC_BO, due to the assumption of partition equilibrium. The standard permeability equation to predict the passive transport across the basolateral membrane is

or

where

(2)

_B is simply the whole basolateral volume accessible to the substrate, since K_BOV_BO is the aqueous equivalent volume of the lipid in the basolateral side. A_B denotes the area of the basolateral outer membrane monolayer capable of permeation. Note that P_BCA_B/_B has the units of a first-order rate constant, s^–1.

We assume that the passive permeability coefficient is symmetric across the basolateral membrane into the cytosol, i.e., P_BC = P_CB. It is also the case that substrates can permeate through cell monolayers between the cells, through the tight junctions, which is called paracellular permeation (Ho et al., 2000). Measured passive permeability coefficients contain both terms, and perhaps others. For all our experiments, radiolabeled C¹⁴-mannitol was used to estimate paracellular permeability coefficients, which were <3% of the passive permeability coefficients measured for amprenavir, quinidine, and loperamide (Tran et al., 2004).

Following the same logic as for Eq. 2, the other mass action equations for passive permeability are written as

where

(3)

As with the basolateral membrane, we assume that passive permeability across the apical membrane to the cytosol is symmetric, i.e., P_AC = P_CA. However, the permeability coefficient across the apical membrane, P_AC, need not equal P_BC. We cannot measure these permeability coefficients directly in the cell monolayer, but we will account for their effect on total transport as shown below.

When V_DC_D(0) is the total mols of drug initially added to the donor side, whether basolateral or apical chamber, then mass balance equations require at all times that

(4)

The passive mass action equations, Eqs. 2 and 3, satisfy mass balance trivially.

In Tran et al. (2004), we showed that loss of substrate, defined as drug leaving the donor chamber and not arriving in the receiver chamber in a timely fashion, could be accounted for simply and quantitatively if loss were first order, at least between the experimental time points. Of the three drugs tested, only loperamide showed loss. The loss was not due to metabolism or degradation/hydrolysis (Tran et al., 2004), but rather due to binding to cells or the experimental apparatus. Most importantly, the loperamide loss was first order. The point of this correction was not to discover the mechanism of loss, but rather to calculate the correct passive permeability coefficients when mass balance is a problem. We can add the substrate loss term, using k_v to denote the first order rate constant of loss:

(5)

Whether drug loss occurs in the aqueous phase and/or in membranes is irrelevant, due to the equilibrium partition assumption.

When loss is first order, we can independently fit the value of the loss rate constant using the average drug concentration over the entire system, including cell plasma membranes. This is defined by

where

(6)

V_T is the effective whole system volume, including partitioning into the cell membranes.

We can now show that the average concentration depends only upon the drug loss rate constant, i.e.,

(7)

using Eqs. 5 and 6. The solution to Eq. 7 is

(8)

Thus, we can also measure total loss of substrate and correct the estimation of passive permeability coefficients for that loss, as done in Tran et al. (2004).

Finally, we need to consider the issue of estimating membrane areas and the fitting of permeability coefficients. In the passive permeation experiment, we fit the product of the permeability coefficient and the area of membrane access. In the preceding equations, A_A and A_B denoted the area of the apical and basolateral membranes capable of permeation. The question is: what area do we use? A morphometric electron microscopy study has suggested that the total basolateral and apical membranes of MDCKII cells grown on polycarbonate filters, as we use here, have roughly the same area, and that this area is 8 times the cross-sectional area of the cell (Butor and Davoust, 1992). However, the actual membrane area of access for permeation is unknown, e.g., due to apical membrane invaginations and basolateral membrane tight binding to the filter support.

It is customary to use the area of the insert onto which the cells grow as the reference for accessible area, i.e., the accessible area would be some multiple of this cross-sectional area. We have assumed that both membranes had areas of access two times that of the cross-sectional area of the insert, i.e., like spheres attached at their equators by the tight junctions. The reason for this choice was simply to align the reported values of the permeability coefficients. Typically, permeability coefficients are fitted assuming a single barrier equation, as was done in Tran et al. (2004). The same data fitted by a single static barrier model of area A or by a two barrier model, each barrier with an area of 2A, will yield the same value for the passive permeability coefficient. Since the bottom line of the passive permeation kinetic analysis is to subtract its contribution to the total transport, this approach is accurate and avoids the confusion of having the same specific permeability represented by two different numbers. Once the issue of membrane access area is resolved, which will not be simple, corrected passive permeability coefficients can be easily predicted from our fitted values. This choice has no effect on the fitting of the active transport parameters.

The second choice we made was how to estimate the values of P_BC and P_AC, given that we could not make intracellular measurements. We fit overall P_BA and P_AB coefficients using the one barrier model (Tran et al., 2004), i.e., the average value across both cell membranes in the B > A or A > B directions, respectively. Depending upon drug, Tran et al. (2004) found that P_BA and P_AB did not have to be time independent or equal. The overall passive permeability barrier of the cells depended upon which side was the donor. To correct the total transport for the contribution by the passive transport, which is the goal in the first place, we have chosen to set P_BC = P_BA and P_AC = P_AB. This allows the donor face to dominate the estimated permeation of drug into the cell, as found in Tran et al. (2004). If the passive permeability coefficients were not constant in time, which was the case for quinidine and loperamide, then we best fit the values for each time interval measured and used an "appropriate" value for each time interval for fitting the active transport parameters, as explained in detail below.

The active transport kinetic model
P-gp is generally assumed to follow the standard Michaelis-Menten reaction (Stein, 1997; Ho et al., 2000; Ambudkar et al., 2003). This reaction takes place within the apical membrane inner monolayer:

(9)

where T₀ is the empty transporter, C_PC the substrate in the apical membrane inner monolayer, T₁ is the transporter bound by substrate, and C_A the substrate after efflux into the apical chamber.

The nomenclature used here, e.g., k₂₁, is read as the second reaction for site 1. The first reaction is association to P-gp and the second reaction is efflux into the apical chamber. This nomenclature was developed to allow for the possibility of two substrate binding sites, such as has been suggested by equilibrium dialysis binding for P-gp (Martin et al., 2000) and for the bacterial multi-drug resistance transporter LmrA in Lactococcus lacti (van Veen et al., 1998). More sophisticated models, including explicit ATP hydrolysis steps, have been proposed (Senior et al., 1995; van Veen et al., 2000; al-Shawi et al., 2003), based upon equilibrium dialysis binding studies, ATP hydrolysis, and/or Michaelis-Menten steady-state kinetic analysis. We shall see for P-gp that the rate constant of association, k₁₁, is so large that whether there is one site or two roughly identical sites, both sites would be filled or empty. Functionally, this means that the number of sites is not kinetically relevant for single substrate measurements. Future experiments examining substrate competition can test the number of binding sites directly. Because of the plan to examine competition, we have retained this notation. In addition, functional two-site kinetics may apply to other multi-drug resistance transporters.

When we append the active transport mass action reaction to the passive permeability kinetics, we obtain the comprehensive kinetic equations for P-gp mediated transport:

(10)

Since we assume that the transporter captures the drug from the apical membrane inner monolayer, the relevant volume is V_AO and the relevant partition coefficient is K_PC. It might appear that the volume for the receptor is 2V_AO, since the transporter crosses the membrane, but its active site is thought to be restricted to the apical membrane inner monolayer (Sharom et al., 2001; al-Shawi et al., 2003), which technically has the volume of V_AO. With respect to fitting parameters, the difference is not significant.

There are several parameters here, but all are important and necessary to obtain realistic estimates for the three elementary rate constants (k₁₁, k_r1, k₂₁) and the density of efflux active P-gp in the apical membrane, denoted T(0). Again, passive permeability coefficients and rate constants for loss are fitted in a separate experiment using the potent P-gp inhibitor GF120918. Partition coefficients are estimated using liposomes whose compositions roughly mimic the cell membrane monolayers.

We can now complete the description of correcting the total transport curve for the effects of passive permeability and loss of drug. Tran et al. (2004) showed for amprenavir that the fitted passive permeability coefficients, P_AB and P_BA, change up to the first half hour and then stabilize. There was no significant loss of amprenavir, but the calculated values for k_v were noisy for the first hour. The other drugs showed different transients for P_AB, P_BA, and k_v. Since the purpose of these fittings was to remove the passive permeability from the total transport, we needed to take these transients into account.

We assume that the passive permeability is the same whether or not P-gp is inhibited by GF120918, so the P_AB and P_BA values are taken in the presence of GF120918. It is nearly impossible to test this assumption, but it is reasonable. During active transport, the true system average drug concentration would have to include the term for T₁, since that is drug-bound to P-gp. When that is done in Eq. 6, then Eq. 7 would contain dT₁/dt and the sum gives the same answer as Eq. 8. Thus, we can estimate the k_v when P-gp is active using Eq. 8. For loperamide, this makes a difference, since Tran et al. (2004) showed that loperamide loss was substantially greater in the presence of GF120918. The reason for this is not known yet.

To incorporate the transients of passive permeability coefficients and k_v values requires a choice. One extreme would be to use all of the individual measurements for P_AB, P_BA, and k_v for each time interval. This would be correct if the passive permeability transients measured in the presence of the P-gp inhibitor GF120918 were exactly the same as in its absence, where we measure total transport. This seems unlikely and we did not want to be fitting the noise of the experiment. The other extreme would be to use the time average of all of the individual measurements for P_AB, P_BA, and k_v. This would largely ignore the transients.

Our approach was to use the individual values until they roughly stabilized and then use the time average from that point on for fitting the remaining time intervals. However, we did not want to make the decision for the stabilization times for each data set, since that would introduce an additional fitting decision. So we used the times found in Tran et al. (2004). For both amprenavir and quinidine, we used individual values for k_v up to 1 h, and the average k_v from all points after 1 h. For the passive permeability coefficients, amprenavir was stable at half an hour, whereas quinidine showed changes for 3 h. For loperamide, both passive permeability coefficients and drug loss coefficients changed for the first 3 h before stabilization.

	ALGORITHM FOR HIERARCHICAL EXHAUSTIVE DATA FITTING

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS KINETIC MODEL OF TRANSPORT... ALGORITHM FOR HIERARCHICAL... RESULTS DISCUSSION APPENDIX: GLOSSARY ACKNOWLEDGEMENTS REFERENCES

 
We now use the fitted passive permeability parameters to find all active transport parameter vectors (T(0), k₁₁, k_r1, k₂₁) that will best fit all of the data for total transport. We start by creating a grid of initial vectors of parameter values that spans all possible values for each of the four active transport parameters. The problem of searching a large multi-variable parameter space for kinetic best fits is there are likely to be very many peaks and valleys in the CV surface, just like the protein folding problem (Chan and Dill, 1998). Search routines like fminsearch stop in the first valley encountered. To illustrate this fitting surface, we use a cartoon with two parameters and the oversimplified fitting surface in Fig. 2, since the true fitting surface for our problem depends upon four parameters and cannot be clearly visualized. The initial parameter vector grid we set up can be represented by the grid on the fitting surface, as projected up from the kinetic parameter (X, Y) plane. For each initial vector, MATLAB numerically integrates the kinetic equations, Eq. 10, computes the CV between the data and the simulation, and then the MATLAB function fminsearch adjusts the parameters until a best local fit to the data is obtained. MATLAB's fminsearch stops when a local minimum in CV is reached within the tolerances of the Nelder-Mead simplex (direct search) method on which this minimization routine is based. The consensus answer we seek will be the cluster of fit vectors that are best fits for all drug concentrations used.

View larger version (51K): [in this window] [in a new window]   FIGURE 2  The general problem of exhaustive multi-parameter fitting is illustrated by a simple theoretical fitting surface over just two fitting parameters, because our model has four fitting parameters and cannot be clearly visualized. Because the fitting surface can have many local minima, a very wide initial net of initial conditions must be cast to assure that the global best fits are not missed. The real shape of this surface depends upon the model and the data.

Initial parameter grid
For the rate constants, we begin with plausible ranges, but these parameters are unconstrained, so that fminsearch in MATLAB can move out of the initial range if there is a local minimum for the CV out there. The rate constants are only constrained to be positive.

The initial range for the association rate constant, k₁₁, was 10⁶–10¹² M^–1s^–1, the upper bound probably being above the lipid lateral diffusion control (Keizer, 1987; Molski et al., 1996; Hinterdorfer et al., 1997). The initial range for the efflux rate constant, k₂₁, was 1–10⁴(s^–1), which contains the range of currently known ATP hydrolysis rate constants (Urbatsch et al., 2003; Loo and Clarke, 2000). The initial range for the dissociation rate constant, k_r1, was actually fixed by the substrate binding to P-gp, given by the binding constant K_B = k₁₁/k_r1, which will start in the typical range of 10²–10⁶ M^–1 . Again, these are just starting ranges and the fitting was not constrained to remain within these ranges.

The amount of P-gp per cell has been estimated for other cell lines using quantitative Western blots (Ambudkar et al., 1997), but not for the cells we use. Since the expression levels of P-gp depends upon many factors, including some that are unknown at this time, we believe that it is important not to fix this parameter. We refer to T(0) as the density for efflux active P-gp. The terminology used is not meant to imply that other P-gp in the cell are misfolded or sequestered within the cell, but our kinetic analysis only measures those P-gp that efflux drug directly into the apical chamber. We will refine this idea below.

On the other hand, the surface density for any membrane protein has physical constraints. It must be less than the surface density of proteins known to closely pack on the membrane surface, e.g., influenza hemagglutinin (HA) on a viral surface, which is a homotrimer with three transmembrane domains and densely covers the viral surface. In comparison, P-gp has 12 transmembrane domains, which clearly constrains P-gp to a smaller surface density than influenza HA. For computational convenience, T(0) was expressed in mols P-gp/L(apical membrane inner monolayer), but this can be converted to P-gp/µm² (of apical membrane) simply by multiplying by 1.2 x 10⁶, which assumes a 2 nm thick acyl chain region in the apical membrane inner monolayer. For example, the surface density of influenza HA on the virion is 16,000/µm² (Mittal et al., 2002), which would mean that T(0) < 16,000/1.2 x 10⁶ M = 0.013 M. We have chosen a minimum surface density for P-gp of four orders of magnitude smaller, which amounts to 1 P-gp/µm², which is very conservative. Thus, the range for T(0) has been fixed as 10^–2 M > T(0) > 10^–6 M. It will turn out that the surface density of efflux active P-gp in these cells was smaller than the estimates of P-gp densities obtained from other cell lines using Western blots. In the Discussion, we will suggest why this might be a perfectly consistent result.

These generous ranges are necessary to make sure that all possible parameter vectors that can fit the data are discovered. The problem then becomes that the ranges are so large that the initial search space is huge and cannot be covered in a brute force way, i.e., starting at one end and going to the other end in some linear fashion. Each of the parameters has a theoretical range of at least four orders of magnitude. On a unitary log scale, each parameter generates nine values per decade, i.e., 1–9, with four decades per parameter, yielding (9 x 4)⁴ = 1.7 x 10⁶ initial vectors. Each initial vector was the starting point for a local best fit, which requires at least 5 min on a 2.4 GHz PC. Brute force would require more than 15 years of computation per drug concentration. So we must build a fitting algorithm.

The hierarchical fitting algorithm
Previous experience in fitting kinetic equations for liposome aggregation and fusion (Bentz and Nir, 1981; Bentz et al., 1983), showed that the fitted ranges for the "forward" parameters were not likely to be very sensitive to the reverse reaction. Since those differential equations were similar to the mass action kinetic equations for the Michaelis-Menten reaction, we took that as the first step of the fitting algorithm. That is, as a start, the "forward" B > A data could be fitted using only the "forward" parameters T(0), k₁₁, and k₂₁, with k_r1 fixed at zero. This reduces the first search to three parameters instead of four, which saves a lot of time since fitting time increases roughly by the power of the number of parameters being fitted. Once the forward parameters were fitted to reduce their ranges, then the "reverse" A > B data could be fitted using k_r1 and the reduced range of values found for T(0), k₁₁, and k₂₁. Of course, if this assumption were wrong, the fit to the A > B data would find very different parameter ranges for T(0), k₁₁, and k₂₁. We shall see that this approach works very well and reduces the numbers of initial grid vectors from more than a million to only a few thousand, requiring <2 weeks of computation per drug concentration. The algorithm also produced an interesting artifact, which we will discuss and discard below.

"Model data" as a control
An extremely important control for fitting several parameters to a kinetic model is to simulate "model data" and subject that "model data" to the same analysis as the real data. "Model data", for the Michaelis-Menten mass action model we have used, was simulated using a parameter vector from the "center of the box" of the consensus vectors for the amprenavir data, as obtained at the end of our data-fitting process. If the transport of amprenavir by P-gp is modeled adequately by the one-binding site Michaelis-Menten reaction, then subjecting the "model data" to the same fitting algorithm must yield the same basic results as observed for the amprenavir data.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS KINETIC MODEL OF TRANSPORT... ALGORITHM FOR HIERARCHICAL... RESULTS DISCUSSION APPENDIX: GLOSSARY ACKNOWLEDGEMENTS REFERENCES

 
The passive and active transport of 90 µM amprenavir across the MDCKII-MDR1 cell monolayers is shown in Fig. 3. The passive permeability was determined in the presence of 2 µM GF120918, a potent inhibitor of P-gp (Hyafil et al., 1993; Polli et al., 2001; Tran et al., 2004), shown by solid symbols. The nmol transported is symmetric over time, i.e., the same for B > A and A > B, which would be expected for a static passive barrier. That was not the case for quinidine and loperamide (Tran et al., 2004). The open symbols show the total transport, when P-gp is active (i.e., without GF120918), with triangles showing the P-gp mediated activation of transport in the B > A direction and the inhibition of transport in the A > B direction. For amprenavir, the passive permeability coefficients were symmetric and constant in time after 15 min, i.e., P_BA = P_AB 200 nm/s. Amprenavir had no significant loss of substrate, i.e., k_v < 1 x 10^–6 s^–1.

Above, we have claimed there are only two kinetic barriers for transport of the drugs studied here, i.e., the basolateral and apical membranes. This means that a third kinetic barrier between the apical and basolateral chambers, e.g., through the cytosol, is kinetically irrelevant for these drugs. We can now prove this assertion.

Relative to passive permeability, the activation of B > A transport by P-gp, i.e., from the basolateral chamber to the apical chamber, illustrated in Fig. 3 by the length of the arrow aimed up, was always greater than the inhibition of A > B transport by P-gp, i.e., from the apical chamber to the basolateral chamber, illustrated by the length of the arrow aimed down. Remarkably, this simple test is diagnostic for there being only two barriers to active transport. If there were an additional kinetic barrier between the basolateral membrane inner monolayer and the apical membrane inner monolayer, e.g., substrates following a path through the cytosol, then the predicted shape of Fig. 3 would be just the opposite. The activation arrow would be shorter than the inhibition arrow (simulations not shown). The reason for this is that the additional barrier slows access of substrate to P-gp from the basolateral side, reducing activation, whereas transport from the apical side is hardly affected, since drug escape from P-gp into the apical membrane is rate limiting anyway. Other drugs have shown this "two-barrier" shape (Troutman and Thakker, 2003). Of course, this diagnostic can identify substrates for which there is a third barrier.

Thus, it appears that the nearly all of B > A transport of amprenavir, as well as loperamide and quinidine, follows the path of binding to outer basolateral monolayer, flip-flopping across the bilayer, and then rapid lateral diffusion within the lipids to the apical membrane inner monolayer, where it binds to P-gp, in a highly reversible way, and is effluxed out into the apical chamber. A > B transport is simply the opposite.

Passive and active transport data were collected for 100, 80, 70, 60, 50, 10, and 1 µM amprenavir, just as was done for 90 µM amprenavir in Fig. 3. We did not use the data from 1 and 10 µM amprenavir for fitting parameters and they became a tester data set at the end. The first fitting step is a coarse fitting the B > A data using only the "forward" parameters, with k_r1 = 0. We set up a grid of initial vectors covering a four-order of magnitude range for each parameter. The grid and tolerance is shown in Run I of Table 1. All possible combinations of the three parameters were used as initial vectors, 729 in all.

View larger version (11K): [in this window] [in a new window]   FIGURE 4  Coefficients of variation between the data, CV, and the 729 Run I fits to the 90 µM amprenavir data are plotted in a rank order, from best, CV 0.03, to worst, where CV > 0.4. Fits with a CV < 0.05 are visually very good. The rank-ordered first 609 fits essentially give the same curve, but from a wide range of parameters. Subsequent rank-ordered fits tracked increasingly worse fits, because MATLAB fminssearch was stuck in a local minimum. Fits to the right of the cut were discarded.

View larger version (33K): [in this window] [in a new window]   FIGURE 5  Fitted pairs of the association rate constant, k11, and P-gp density, T(0), for 90 µM amprenavir remaining after the cut in Fig. 4 are graphed here. The range of the parameters has not been reduced in the fitting and they are uncorrelated, except for a small trend line at the bottom, which is examined in Fig. 8.

View larger version (33K): [in this window] [in a new window]   FIGURE 8  Final forward fits for the amprenavir data. Panel A shows all pairs of T(0) and k11 that best fit the final fitting of the forward B > A amprenavir data, Run III. The initial grid for T(0) and k11 is the box with a dashed line border. The amprenavir concentrations are color coded, as shown by the legend. The box labeled "consensus fits" contains pairs from all amprenavir concentrations, except 50 µM, which had only correlated fits. Thus, any pair within the box gives essentially the same best fit to all amprenavir concentrations. There is also a series of fits labeled "correlated fits", the straight lines at the bottom. These fits are different for each amprenavir concentration and show no reduction in volume. As explained in the text, these fits are artifacts from applying the fitting algorithm to the mass action kinetics of the Michaelis-Menten reaction. These fits may be discarded. Panel B is a blowup of the consensus fits, showing the distribution of fitted pairs. We shall see in Fig. 11 below that using parameter values in the center of this box gives good fits to all amprenavir concentrations, despite the fact that 1, 10, and 50 µM amprenavir yielded only correlated fits (data not shown).

View larger version (16K): [in this window] [in a new window]   FIGURE 6  Fitted pairs of the efflux rate constant, k21, and the P-gp density, T(0), for 90 µM amprenavir remaining after the cut in Fig. 4 are graphed here. The possible range of these parameters has not been reduced in the fitting, but there is essentially complete correlation. Thus, for 90 µM amprenavir, the product of k21T(0) 0.01 M/s.

To show that this correlation is not an artifact of the fitting algorithm, we tested it with the "model data". We generated computer-simulated data using a parameter vector within the set of best fits for amprenavir, as defined in Table 4 below. These "model data" represent perfect data for the case that P-gp is a single site Michaelis-Menten type transporter and there is no experimental error. We note here that the "model data" for P-gp with one or two binding sites were the same, assuming both sites had to be filled for one drug molecule to be effluxed. This is because association was so rapid (simulations not shown). The "model data" for all drug concentrations were then fitted with the algorithm being followed here. For each concentration of substrate, the fitting gave a constant value for k₂₁T(0), which is shown as the solid hyperbolic line in Fig. 7. Interestingly, the small substrate concentrations predict a k₂₁T(0) product that is smaller than the value of k₂₁T(0) used to simulate the data in the first place, i.e., 9 x 10^–3 M/s. However, as the substrate concentration increased, the predicted k₂₁T(0) becomes equal to the correct value. Thus, even with model data, the correct k₂₁T0 value is obtained only when the substrate concentration is large enough. This fit was the outcome of Run I in Table 1.

View larger version (15K): [in this window] [in a new window]   FIGURE 7  Fitted values of k21T(0) for 50–100 µM amprenavir, (•), and predicted from simulated "model data" (solid line). Model data were simulated for the Michaelis-Menten reaction model for P-gp, without experimental error, using the parameters given in Table 4 for amprenavir. At each concentration, a predicted k21T(0) is found, just like in Fig. 6. The solid hyperbolic line is the predicted k21T(0) for all the model data. At high substrate concentrations, the k21T(0) used to simulate the model data was recovered. However, at low concentrations, the predicted k21T(0) was smaller than the correct value. The k21T(0) values predicted for each concentration of the amprenavir data are shown by solid circles. For each concentration, the best 300 fits gave the same value for k21T(0), i.e., the standard deviation for each point is much smaller than the size of the symbols. The average of the higher amprenavir concentrations was roughly k21T(0) = 0.009 M/s, which is the value shown in Tables 2 and 4 as k21T(0) = 180 (1/s) x 50 (µM) = 0.009 M/s. Thus, for all subsequent fittings for amprenavir, for each value of T(0) tested, we fixed the value of k21= (0.009 M/s)/T(0).

The second or fine fitting comes next. As explained in the Fitting Algorithm section, this is a two-parameter fit to the B > A data using T(0) and k₁₁, with the parameter range of the first, or coarse fit, but with a finer grid. Also, k_r1 = 0 and k₂₁ = 0.009/T(0) s^–1. A higher stringency is used with ode23s; relative and absolute tolerances were set at 10^–9. The grid and tolerance is shown in Run II of Table 1. All possible combinations of the two parameters were used as initial vectors, 13² = 169 in all. For economy, these fits are not shown. The parameter ranges were reduced to 3