Design of Optimized Hypoxia-Activated Prodrugs Using Pharmacokinetic/Pharmacodynamic Modeling

Hypoxia contributes to resistance of tumors to some cytotoxic drugs and to radiotherapy, but can in principle be exploited with hypoxia-activated prodrugs (HAP). HAP in clinical development fall into two broad groups. Class I HAP (like the benzotriazine N-oxides tirapazamine and SN30000), are activated under relatively mild hypoxia. In contrast, Class II HAP (such as the nitro compounds PR-104A or TH-302) are maximally activated only under extreme hypoxia, but their active metabolites (effectors) diffuse to cells at intermediate O2 and thus also eliminate moderately hypoxic cells. Here, we use a spatially resolved pharmacokinetic/pharmacodynamic (SR-PK/PD) model to compare these two strategies and to identify the features required in an optimal Class II HAP. The model uses a Green’s function approach to calculate spatial and longitudinal gradients of O2, prodrug, and effector concentrations, and resulting killing in a digitized 3D tumor microregion to estimate activity as monotherapy and in combination with radiotherapy. An analogous model for a normal tissue with mild hypoxia and short intervessel distances (based on a cremaster muscle microvessel network) was used to estimate tumor selectivity of cell killing. This showed that Class II HAP offer advantages over Class I including higher tumor selectivity and greater freedom to vary prodrug diffusibility and rate of metabolic activation. The model suggests that the largest gains in class II HAP antitumor activity could be realized by optimizing effector stability and prodrug activation rates. We also use the model to show that diffusion of effector into blood vessels is unlikely to materially increase systemic exposure for realistic tumor burdens and effector clearances. However, we show that the tumor selectivity achievable by hypoxia-dependent prodrug activation alone is limited if dose-limiting normal tissues are even mildly hypoxic.

In the present study we utilize pharmacokinetic/ pharmacodynamic (PK/PD) models to explore relationships between the reaction/diffusion properties of HAP in the tumor microenvironment and their antitumor activity and selectivity. In this context it is useful to distinguish two broad classes of HAP with different PK/PD features (Figure 1). Class I HAP are activated by reduction under relatively mild hypoxia to generate a reactive cytotoxin, which is restricted to the cell in which it is formed. This class is typified by the benzotriazine N -oxides tirapazamine and SN30000 that undergo 1-electron reduction to short-lived cytotoxic radicals (27,28), which appear not to affect adjoining cells in three-dimensional (3D) cell cultures (29). The O 2 concentration for 50% inhibition of cytotoxicity (K O 2 ) in stirred single cell suspensions is approximately 1 µM for both tirapazamine (30,31) and SN30000 (32). Class II HAP (typified by nitro compounds such as PR-104A and TH-302) require more severe hypoxia for 1electron reduction, but the resulting radicals are further reduced to a relatively stable effector that can diffuse from the cell of origin (29,(33)(34)(35); this diffusion of effector to nearby cells (which may not themselves be capable of prodrug activation) is here referred to as a bystander effect. In the case of PR-104A, a K O 2 value of 0.1 µM has been estimated in stirred single cell suspensions (31), and activation of TH-302 also appears to require severe hypoxia (34) as is the case for other nitro compounds (36)(37)(38). Given that half-maximal radiosensitization of tumor cells by oxygen requires concentrations of~4 µM (39,40), the ability of Class II HAP to overcome radioresistance may depend on bystander effects. We have suggested (41) that Class II HAP may be preferable to Class I because activation will be confined to the extreme hypoxia found in tumors, thus minimizing toxicity to normal tissues with mild, physiological hypoxia [e.g., retina (42), liver (43)(44)(45), esophagus (46), skin (47,48), and possibly the bone marrow stem cell niche (49) although the oxygenation status of the latter is controversial (50)]. The bystander effect from Class II HAP may contribute to the reported monotherapy activity of PR-104 (33,51,52) and TH-302 (53) in preclinical models. However, it is not known under what conditions this theoretical advantage for Class II HAP might be realized, or how this activity might be optimized by prodrug design, or whether diffusion of active metabolites into the tumor microvasculature might contribute to systemic toxicity if this process is too efficient.
Spatially resolved pharmacokinetic/pharmacodynamic (SR-PK/PD) models provide tools for addressing these questions. These models can be used to describe concentration gradients of oxygen, HAP, and their effectors in tumors, using mapped microvascular networks, and to calculate resulting reproductive cell death (clonogenic cell killing). These models include the effects of heterogeneity in inter-capillary distances, vessel diameters, blood flow rates, and vessel oxygen and drug concentrations. We have validated an SR-PK/PD model for Class I HAP by showing that it predicts activity of tirapazamine analogs combined with radiotherapy in human tumor xenografts (32,54,55). This modeling clearly demonstrated the need to optimize rates of reductive metabolism such that penetration into hypoxic regions is not compromised by excessive consumption of the prodrug. Recently, we have also reported an SR-PK/PD model for the Class II HAP, PR-104, and used this to estimate that 30-50% of its activity in FIGURE 2 | Schematic representation of a generalized HAP PK/PD model. Transfer of prodrug and effector between the extracellular and the intracellular compartment is defined by rate constant k t N , where N refers to each compound. In the extracellular compartment compounds can diffuse as defined by their diffusion coefficients D N (double-headed arrows). Prodrug activation (defined by k met P ) is restricted to the intracellular compartment and is O 2 -dependent (see Eq. 3). For an effector that confers cytotoxicity by way of reacting irreversibly with its target (case 1, e.g., DNA-alkylating agents) the model assumes that potency a scales with reactivity of the effector (k met E ), while in case 2 (e.g., a reversible inhibitor) a is independent of k met E .
HCT116 and SiHa xenografts is due to bystander effects both as monotherapy and combined with radiation (35). Here, we use a generalized SR-PK/PD model, in which a HAP is metabolized by an oxygen-inhibited process to a single effector (Figure 2), to ask under what conditions Class II HAP might provide greater tumor activity and selectivity than Class I HAP, and to identify the prodrug features required for optimal antitumor activity. This generalized HAP model makes explicit the diffusion of both the prodrug and effector in the extracellular (interstitial) compartment. We consider two types of effector, which elicit cytotoxicity via irreversible reaction with a target (Case 1, as for an alkylating agent) or by reversible binding to its target (Case 2, as for a receptor ligand).

METHODS
The generalized SR-PK/PD model calculates steady-state concentrations of oxygen, HAP, and effector as well as resulting cell killing in digitized 3D tissue microregions using Green's function methods (35,54,56). We used two different tissue microregions that were derived by mapping microvascular anatomy as well as direction and velocity of blood flow in a rat cremaster muscle (56) ("normal" network) and a subcutaneous FaDu tumor xenograft (57) ("tumor" network). The blood vessels are represented by cylindrical segments and vessel walls are treated as part of the tissue space. The model was implemented using a customized version of the Green's function method written in Visual C++ (Microsoft Visual Studio 2010 Express) (35,58).

CALCULATION OF OXYGENATION
Convective transport of oxygen along vessel segments and diffusion into the surrounding tissue (represented as homogeneous Frontiers in Oncology | Pharmacology of Anti-Cancer Drugs medium) was calculated based on previous estimates for blood flow in the networks, blood O 2 content, tissue diffusion, and consumption of O 2 (35,56). All O 2 transport parameters are summarized in Table 1.

CALCULATION OF PHARMACOKINETICS
Inflow of prodrug to the microvascular networks was defined by its unbound plasma level, using area under the concentration-time curve (AUC) as a time-independent exposure variable compatible with Green's function formalism. Vessel walls were modeled as offering negligible resistance to radial flux by setting the intravascular resistance constant K (56) to a low value (0.1 s/µm). Extravascular transport was calculated using a (pro)drug transport model (Figure 2) assuming that the tissue consists of two homogeneous (extracellular and intracellular) compartments, with activation of prodrug restricted to the intracellular compartment and diffusion confined to the extracellular compartment. Concentrations in the two compartments were calculated as follows: C e and C i are the extracellular and intracellular concentrations, respectively, of prodrug or effector (denoted by N = P or E), ϕ i and ϕ e are the intra-and extracellular volume fractions with ϕ e = 1 − ϕ i , k tN are first order rate constants for cell transfer between the extracellular and intracellular compartments, D N is the effective diffusion coefficient of compound N, 2 is the Laplacian operator, k metN is the rate constant for drug metabolism and r is the rate of metabolic production with r P = 0 and rate of effector production equal to the rate of loss of prodrug: r E = ϕ i k metP C iP . The rate constant for prodrug activation is O 2 -dependent and is given by: where k metP,max and k metP are the prodrug activation rate constants under anoxia and at the O 2 concentration [O 2 ], respectively, and K O 2 is the O 2 concentration for half-maximum prodrug activation.

CALCULATION OF CELL KILLING
Survival probability (SP) at each point of the tumor microregion was calculated using the PK/PD relationship: where AUC iE is the intracellular exposure (area under the intracellular concentration-time curve) to the effector and a is a proportionality constant (potency) equal to the reciprocal of the effector AUC to produce a SP of 0.1 as measured by clonogenic cell survival. We distinguish two cases: case 1, where potency a scales with k metE , broadly representing irreversible inhibitors such as DNAalkylating agents for which an increase in reactivity has similar effects on reaction with target and non-target nucleophiles (59,60) and case 2, where a is independent of k metE (broadly representing reversible inhibitors). Case 1, where cell killing is proportional to the reactivity of the effector is used throughout this study unless noted otherwise.
To estimate cell SP after radiation treatment, the linearquadratic model was used as previously (35,54). The SP from  (56).
These O 2 transport parameters gave hypoxic fractions and median pO 2 values (Figure 3) in the range measured for tumors and normal tissues in humans using O 2 electrode histography (45).
www.frontiersin.org both radiation and prodrug at each point of the tumor microregion was calculated from the sum of log (SP) due to drug and radiation alone. Log cell kill was then calculated as log cell kill = − log(Average SP) (5) where the SP is averaged over the whole tumor microregion. Prodrug-induced cell kill in addition to radiation was calculated as the difference between overall log cell kill by prodrug + radiation and log cell kill by radiation alone.

ESTIMATION OF NET TRANSPORT ACROSS THE TUMOR-PLASMA BOUNDARY
The SR-PK/PD model calculates the transport of prodrug and effector across the walls of the 582 vessel segments in the FaDu tumor microvascular network. Summing the 582 values for each compound gives the net transport across the tumor-plasma boundary (in picomoles). This was normalized to the volume of the FaDu tumor microregion.

PHARMACOKINETIC STUDIES IN MICE
All animal studies were approved by the University of Auckland Animal Ethics Committee. Male NIH-III nude mice (NIH-Lyst bg Foxn1 nu Btk xid ; 18-20 g body weight), derived from breeding mice purchased from Charles River Laboratories (Wilmington, MA, USA), were dosed i.v. with PR-104H. The dosing solution was prepared fresh for each mouse by dilution of a stock solution of PR-104H in acetonitrile with sterile saline (to <20% acetonitrile) within 3 min of dosing to minimize loss of PR-104H in aqueous solution. At the dose volume of 0.005 ml/g the acetonitrile dose was 755 mg/kg and the vehicle alone did not cause signs of toxicity over the duration of the experiment. Blood was collected by cardiac puncture up to 3 min after cervical dislocation as in (62). Plasma was prepared as described (63) and mixed with three volumes of cold acidified methanol (methanol:ammonium acetate:acetic acid 1000:3.5:0.2, v/w/v)~6 min after cardiac puncture. Liver tissue was dissected within 3 min of blood sampling, snap-frozen in liquid nitrogen, and stored at −80°C along with plasma samples. Frozen tissue was pulverized at liquid nitrogen temperature using a BioPulverizer™ (BioSpec Products, USA), transferred to tared tubes, weighed, vortex mixed for 30 s with an equal volume of cold phosphate buffered saline and six volumes of cold acidified methanol, and stored at −80°C. For LC-MS/MS analysis [as reported in Ref. (35)], tissue and plasma extracts were centrifuged (13,000 × g, 4°C, 10 min) and supernatants were mixed 1:1 with cold water containing 1 µM PR-104H-d4 and 1 µM PR-104M-d4.

RESULTS
In order to identify optimal properties of HAP, we developed a SR-PK/PD model that captures the key features of both Class I and Class II HAP (Figure 2). This model is a generalized form of our SR-PK/PD model for PR-104 (35), which, unlike earlier SR-PK/PD models (31,32,54) accommodates diffusion and reaction of the effector and explicitly considers intra-and extracellular compartments in order to represent the plasma membrane barrier to prodrug uptake and effector efflux. Our generalized model assumes that the effector can be further converted to inactive products within the cell (defined by rate constant k metE ). The model incorporates only one effector (unlike the two alkylating metabolites from PR-104A) and represents an idealized HAP that is not a substrate for O 2 -insensitive two-electron-reduction, which is an off-target activation mechanism of some HAP including PR-104A (51). We also ignore any systemic generation of effector outside the modeled tissue region, such as from hepatic metabolism.

TISSUE MICROREGIONS USED FOR SR-PK/PD MODELING
We used a digitized 3D tumor microregion ( Figure 3A) that was derived by mapping a region of a subcutaneous FaDu tumor xenograft grown in a mouse dorsal window chamber (57). To assess HAP activation under conditions of physiological hypoxia we also used a 3D region of a rat cremaster muscle (56) as an example of a normal tissue with a well-organized and efficient microvascular network ( Figure 3B). While the maximum distance to the nearest blood vessel is 50 µm in the cremaster muscle region, 20% of tissue points in the FaDu tumor region are situated further from blood vessels (50-100 µm; Figure 3C). The FaDu tumor region also showed a higher heterogeneity in vessel diameters ( Figure 3D), and a lower median blood flow rate ( Figure 3E) at the values chosen for total blood inflow into the regions (Q; Table 1). Q per tissue volume was 2.5-fold lower in the FaDu tumor than in the cremaster muscle microregion, which is realistic for a low-perfused hypoxic tumor area. The value for the tumor microregion falls within the range of mean blood flow values measured in human tumors, which show high variability with values that can be lower, similar, or higher than in the tissue of origin (64,65). A radiobiological hypoxic fraction (<4 µM O 2 ) of 43% was achieved in the FaDu tumor microregion ( Figure 3F) by using O 2 transport parameters previously chosen to model the hypoxic fraction of tumor xenografts (35) ( Table 1). The O 2 transport parameters for the cremaster muscle region were as reported (56), with the oxygen consumption rate adjusted to achieve a radiobiological hypoxic fraction of 7.5% ( Figure 3F) that is consistent with moderate hypoxia as measured in some normal human tissues using O 2 electrode histography (45). However, the muscle region lacks the severely hypoxic cells (defined as <0.13 µM O 2 , i.e., below the K O 2 value for PR-104A) that comprise 5.0% of the FaDu microregion ( Figure 3F).

MODULATION OF KINETICS OF PRODRUG ACTIVATION
Our previous SR-PK/PD model for Class I HAP showed that high prodrug activation rates limit killing of hypoxic cells by impairing tissue penetration of the prodrug (54). To investigate whether this also applies to HAP with bystander effects and/or activation restricted to more severe hypoxia (low K O 2 ), we modulated the Frontiers in Oncology | Pharmacology of Anti-Cancer Drugs prodrug activation rate constant for anoxia (k metP,max ) under these conditions. For simulations without a bystander effect, k tE was set to 0, which traps effector within the intracellular compartment. Increasing k metP,max resulted in a decrease of prodrug penetration to the most hypoxic regions, an effect that was less pronounced with a low K O 2 ( Figure 4A) than with a high K O 2 value (Figure 4D), as expected. Accordingly, killing in the most hypoxic regions (<0.13 µM O 2 ) was higher and less affected by an increase in k metP,max for low-K O 2 HAP (Figures 4B,E; without bystander effects). The observed shift of cell killing to higher O 2 concentrations with increasing k metP,max reflects increased prodrug activation at intermediate O 2 concentrations. In case of low-K O 2 HAP this improved complementation of radiation-induced killing ( Figure 4B), but in case of high-K O 2 HAP it resulted in a loss of hypoxic selectivity at a 100-fold higher k metP,max ( Figure 4E). If the effector was allowed to diffuse through the plasma membrane (i.e., the bystander model), killing was increased across the whole tumor microregion, and killing <0.13 µM was less affected by decreasing prodrug penetration with increasing k metP,max (Figures 4C,F).
Figures 4C,E illustrate the dependence on prodrug activation rate constant for HAP of Class II (low K O 2 + bystander) and Class I (high K O 2 , no bystander) respectively. It should be noted that in the no-bystander simulations the effector potency, a, was set to ã fivefold lower value (0.0414 µM −1 ) than for Class II in Table 2.
www.frontiersin.org   Figure 4. In all cases, monotherapy activity increased with increasing prodrug activation rates over the range examined (k metP,max = 0.001-1 s −1 ), although more markedly for low-K O 2 compounds (shown in red). This difference was even more pronounced for killing of radiobiologically hypoxic cells (i.e., killing additional to radiation; Figure 5B), in which case activity of high-K O 2 compounds (shown in green) decreases at high rates of metabolism reflecting a severe penetration problem (which is not fully offset by allowing bystander diffusion). Class II HAP (dark red) are thus more tolerant of high values of k metP,max than Class I HAP (light green).
The above analysis suggests that monotherapy activity increases monotonically with k metP,max up to very high values, but it is important to consider whether this changes selectivity relative to normal tissues. To address this we used the above cremaster muscle microregion (Figure 3), to represent a generic wellperfused normal tissue with mild physiological hypoxia. For the present purpose, we assumed that this "normal tissue" is populated by cells with the same intrinsic sensitivity to effector as cells in the tumor microregion. Cell kill in the normal tissue region showed an even steeper dependence on k metP,max than the tumor region ( Figure 5C). As a consequence, tumor selectivity (assessed using the tumor:normal tissue ratio of log cell kill) fell as rates of bioreductive metabolism increased ( Figure 5D). This falling therapeutic ratio reflects the greater compromise to delivery of rapidly metabolized prodrugs in the tumor network because of lower blood flow rates, longer diffusion distances, and more severe hypoxia than in the normal tissue (Figure 3). The model confirmed our expectation that toxicity to physiologically hypoxic normal tissues will be higher for a high-K O 2 HAP (Figure 5C).

MODULATION OF PRODRUG DIFFUSION PROPERTIES
We next asked whether the activity of class II HAP is sensitive to the diffusibility of the prodrug, as previously shown by the SR-PK/PD model for tirapazamine analogs (54). In the latter study, tissue diffusion was defined by a single parameter, the effective diffusion coefficient as measured by flux through multicellular layer cultures; this represents the weighted averages of the free drug diffusion coefficient within cells, across membranes and in the extracellular matrix, based on treating tissue as homogenous space. In contrast, the present model explicitly considers extra-and intracellular compartments, and overall diffusion is determined by two parameters: the rate constant for transfer between extraand intracellular compartments, k tP , and the diffusion coefficient in the extracellular compartment D P (Figure 2). Since changes in the physicochemical properties of HAP are expected to affect membrane permeability more than paracellular diffusion, k tP was modulated. The results (Figure 6) show that hypoxic cell killing for Class I HAP is optimal at k metP,max~0 .1 s −1 [ Figure 6A; consistent with the previous SR-PK/PD model for tirapazamine (54)]. Class Frontiers in Oncology | Pharmacology of Anti-Cancer Drugs I HAP also showed an optimum value of k tP , reflecting the fact that the membrane transfer rate constant controls how much prodrug is available for intracellular activation. In contrast, for Class II HAP killing additional to radiation increased monotonically over the full range of k metP,max and k tP (Figure 6B) reflecting the tolerance of high rates of intracellular prodrug activation owing to a low K O 2 and because effector diffusion offsets compromised prodrug penetration.

MODULATION OF THE EFFECTOR DIFFUSION RANGE
Next, we investigated the dependence of antitumor activity on effector parameters for Class II HAP. Since tissue transport of the effector is dependent on its membrane transfer (k tE ), paracellular diffusion (D E ) and its stability (k metE ), these parameters were modulated~10-fold relative to the base model. In order to see an impact of a change in the paracellular diffusion coefficient D E , k tE had to be increased 10-fold to 0.1s −1 , indicating that effector transport is otherwise limited by membrane transfer. An increase in diffusion coefficient, D E (Figure 7A) caused increased diffusion of the effector away from the hypoxic regions where it is formed, which decreased killing in hypoxic regions and increased killing in aerobic regions. As a consequence, overall killing in addition to radiation decreased and monotherapy activity increased (Figures 7B,C). A similar effect could be observed when increasing the membrane transfer rate constants k tE (Figures 7D,E). However, the effect of increased diffusibility (higher D E or k tE ) on single-agent activity was only minor because the effector was defined to freely permeate the blood vessel walls, and thus lost into the vasculature. The protective effect of this washout was confirmed by simulations with zero vessel permeability to the effector, which showed a much larger increase in HAP monotherapy activity with increasing k tE (Figures 7G,H). Notably, an increase in effector diffusibility in the presence of vessel permeability elevated the tumor:normal tissue ratio of killing (Figures 7C,F), reflecting that washout of effector into the circulation is more protective in normal tissue (with lower distances to nearest vessel; Figure 3) than in tumor tissue. The effect of increasing the rate constant for effector reaction k metE was investigated using the default case where cell killing scales with effector reactivity (case 1) and alternatively where cell killing is independent of reactivity (case 2). The proportionality constant a was set to the values given in Table 2 to achieve equivalent killing of the two effector types at the default k metE of 0.01 s −1 . In case 1, an increase in k metE produces a proportional increase in the rate of formation of cytotoxic lesions and resulted in higher killing in hypoxic regions ( Figure 8A). Killing in aerobic regions (>4 µM O 2 ) was decreased at high values of k metE due to associated low effector penetration, and this limited overall activity ( Figure 8B). Tumor selectivity decreased with increasing k metE ( Figure 8C) due largely to increasing cell killing in normal tissue.
In case 2 (where cytotoxic lesion production is independent of effector reactivity) decreasing k metE resulted in higher killing www.frontiersin.org across the whole tumor microregion (Figure 8D), with a large increase in HAP activity, both as a single agent and in addition to radiation, with increasing effector stability ( Figure 8E). Killing in the normal tissue microregion was affected to a similar extent, so that tumor selectivity was comparable for different values of k metE ( Figure 8F).

ESTIMATION OF THE IMPACT OF EFFECTOR WASHOUT BY BLOOD FLOW
The above analysis demonstrated that washout of effector by blood flow can, under some conditions, significantly protect perivascular cells from cytotoxicity. This raises the question whether effector washout might have pharmacodynamic effects in downstream tumor regions (redistribution of effector through the tumor vasculature), or might result in significant systemic exposure to the effector. To test the sensitivity of effector washout on different parameters, we used the Class II HAP model to calculate net transport of effector across the tumor/plasma boundary in each vessel segment of the tumor microregion. This showed high sensitivity to the rate constant of effector production (k metP,max ) and effector stability (k metE ) but only moderate sensitivity to effector diffusibility (k tE and D E ) (Figure 9), because the former parameters substantially affect overall effector concentrations.
The impact on systemic exposure will depend on the distribution volume (V D ) and clearance (CL) of the effector. To evaluate a specific case, we returned to the SR-PK/PD model for PR-104 (35) because the plasma AUC of its active metabolites (PR-104H and PR-104M) after dosing PR-104 has been determined in three strains of mice (35,52,61). To facilitate this analysis, we measured the pharmacokinetics of PR-104H and PR-104M in plasma and liver following i.v. dosing of NIH-III nude mice with synthetic PR-104H (Figure 10). This showed a very short half-life (2.6 min) for PR-104H in plasma, and high PR-104M concentrations in plasma and liver, indicating rapid conversion of PR-104H to PR-104M. The half-lives of PR-104M in plasma and liver (~7 min) were similar to the half-life of PR-104M in culture medium at 37°C [6.6 min (35)], suggesting that CL of PR-104M is mainly due to its chemical instability.
Since washout was most sensitive to k metP,max (Figure 9), we performed PR-104 SR-PK/PD model simulations using two different values of this parameter, corresponding to rates measured in vitro for HCT116/WT and a cell line with 20-fold higher k metP,max (HCT116/sPOR#6) (35). This confirmed that washout of metabolites produced in hypoxic regions can increase plasma levels in nearby vessels, an effect that was more pronounced at the higher k metP,max (white arrows in Figure 11). Using V D and CL estimates from measured plasma PK following administration of PR-104 or PR-104H (Table 3), we calculated that washout would elevate the free plasma AUC of PR-104H+M in Frontiers in Oncology | Pharmacology of Anti-Cancer Drugs NIH-III nude mice [estimated from the measured PK after i.p. administration of 562 µmol/kg PR-104 (35)] by only 0.8% in case of HCT116/WT tumors, and 3.3% in case of HCT116/sPOR#6 tumors ( Table 3). This was calculated assuming all PR-104H+M was released at once and hence represents a maximum addition to the circulating metabolites. This is consistent with published data showing similar plasma concentrations of reduced metabolites in NIH-III nude mice with HCT116/WT and HCT116/sPOR#6 tumors 30 min after i.p. administration of 562 µmol/kg PR-104 (35).
It should be noted that all simulations in this study were performed at blood flows where prodrug delivery and effector extraction were not highly perfusion limited. To demonstrate this (and the effect of blood flow on increasing effector washout into the systemic circulation) we increased the blood flow rate, Q, 2.5-fold, while decreasing inflow pO 2 , to maintain a similar hypoxic fraction ( Figures S1A,B in Supplementary Material). Under conditions of high prodrug metabolism (k metP,max of 1 s −1 ), this increased prodrug extraction from the plasma 1.5-fold and effector washout into the plasma twofold ( Figure S1C in Supplementary Material) but overall killing in the microregion increased by only 26% (Figures S1D,E in Supplementary Material) reflecting little increase in perivascular effector concentrations.

DISCUSSION
In this study, SR-PK/PD modeling was utilized to identify strategies for optimization of HAP. The model has a solid biological foundation because it is based on our previous SR-PK/PD model for PR-104 for which parameters have been determined experimentally using three cell lines (35). The latter model was the first to explicitly consider the intra-tumor distribution of bystander metabolites.

NORMAL TISSUE TOXICITY
A novel feature of the present study is the comparison of predictions for a tumor and normal tissue network, the latter based on a mapped microvascular network in a rat cremaster muscle (56) with oxygen consumption rate adjusted to simulate mild hypoxia (Figure 3). This comparison is important because the utility of HAP is likely to be constrained by the requirement that prodrug activation at intermediate O 2 does not cause toxicity in tissues with physiological hypoxia, such as liver (43)(44)(45), retina (66), and possibly bone marrow (49). The low tumor:normal tissue ratios of killing determined in our study (Figures 5, 7, and 8) point to limitations on achieving tumor selectivity with HAP. However, it is important to note that our model makes the assumption that intrinsic cellular sensitivity to the released effector (potency parameter a) is the same for cells in the tumor and normal tissue networks. Most effectors of HAP in clinical development are DNAreactive cytotoxins (13) that provide some tumor-specificity due to the cytokinetics and DNA repair abnormalities (67) of tumor cells relative to normal cells. In addition, release of HAP effectors that target oncogenic pathways in tumor cells has potential for conferring additional tumor selectivity (68). Normal tissue models will need to be revisited when more quantitative information about oxygenation and drug targets in potentially dose-limiting normal tissues is available. Nevertheless, the present results emphasize the importance of using HAP strategies to augment selectivity of effectors that already provide some level of tumor selectivity, rather than as the sole mechanism of tumor targeting.
In addition, normal tissue toxicity due to washout of effectors has been a concern in the development of targeted anticancer prodrugs (69,70) although to our knowledge no clear link between toxicity and increased effector levels has been established. In the current study systemic toxicity of HAP as a result of release of active metabolites from the tumor appears unlikely; our model suggests that a significant contribution to toxicity would require a combination of extreme assumptions (high tumor burden, intratumor activation rates, effector stability, tumor blood flow rates, as well as low systemic CL of effector). However, toxicity in normal tissue sharing a microcirculatory system with the tumor remains a possibility.

OPTIMAL PRODRUG PARAMETERS OF CLASS I AND CLASS II HAP
This study compared an idealized Class I HAP (high K O 2 with nobystander effect, illustrated by tirapazamine and SN30000) and a Class II HAP (low K O 2 with bystander effect, e.g., PR-104A and TH-302) to assess which of these HAP strategies is better suited for complementation of radiotherapy. Our previous tirapazamine SR-PK/PD model showed that Class I HAP exhibit an optimum k metP,max , above which hypoxic cell killing decreases www.frontiersin.org FIGURE 7 | Modulation of effector diffusion properties for a class II HAP. SR-PK/PD model simulations using the base model ( Table 2), with modulation of the following effector parameters: (A-C) the extracellular diffusion coefficient D E , using a 10-fold higher value for the membrane transfer rate constant k t E ; (D-I) membrane transfer rate constant k t E , with (D-F) or without (G-I) blood vessel permeability to the effector. Left panels show survival probability (as a function of O 2 in the FaDu tumor microregion) after 15 Gy radiation (gray) or HAP (black: base model; red and blue: 10-fold higher and lower parameter value, respectively). Middle panels show average monotherapy activity (black, solid lines) and killing additional to radiation (white, dashed lines) in the tumor region. Right panels show the tumor:normal tissue ratio of cell killing.
because increasing consumption of the prodrug limits its penetration to hypoxic regions (54). This finding could be replicated using the present model, in which intra-and extracellular compartments are made explicit (Figure 5). In contrast, activity of Class II HAP was predicted to increase monotonically with increasing k metP,max in the investigated range of 0.001-1 s −1 ; we consider this range to be physiologically and pharmacologically relevant based on measured rate constants (32,35,54,60). In addition, the potential to further increase k metP,max by drug design would be limited by membrane transport and by the concurrent decrease in tumor selectivity ( Figure 5D).
The superiority of Class II over Class I HAP at high k metP,max can be explained by two factors. Firstly, the lower K O 2 value minimizes metabolic loss during diffusion into hypoxic target regions (compare Figures 4A,D). Secondly, bystander effects alleviate the impact of poor prodrug penetration due to redistribution of effector between regions of different prodrug exposure. Both of these factors also explain the lower sensitivity of antitumor activity to prodrug diffusibility of class II relative to class I HAP (Figure 6).

DESIGN OF OPTIMIZED CLASS II HAP
Spatially resolved pharmacokinetic/pharmacodynamic modeling represents a valuable tool for the rational design of improved HAP, as has been demonstrated by our previous SR-PK/PD model for Class I HAP (54). The insight that there is an optimum for combinations of D P and k metP that can accommodate the competing requirements of efficient prodrug activation and tissue penetration prompted the SR-PK/PD model-guided screening of tirapazamine Frontiers in Oncology | Pharmacology of Anti-Cancer Drugs  analogs (71) that identified SN30000 as a prodrug with superior tissue penetration and antitumor activity in combination with radiation (32). The model-based design of class II HAP presents a greater challenge as it requires the specific consideration of the transport parameters of the effector in addition to those of the prodrug. In this study SR-PK/PD modeling was used to identify reaction-diffusion parameters of Class II HAPs that could be optimized. The outcome is summarized in Table 4, which shows that antitumor activity of class II HAP could be enhanced by increasing k metP,max , although this would need to be balanced against the need to maintain sufficient tumor selectivity.
When looking at the impact of a change in effector metabolism (k metE ) we distinguished two types of effectors: an effector that needs to react irreversibly with its target to elicit cytotoxicity, meaning that cell killing is a function of this chemical reactivity as for alkylating agents (Case 1), and an inhibitor where chemical transformation of the effector (whether spontaneously or by metabolism) gives rise to non-toxic products (Case 2). Modeling showed that effector reactivity should be increased in Case I but decreased in Case 2, in order to improve antitumor activity. In Case 1 the increase in overall killing when increasing reactivity by drug design needs to be balanced against a predicted decrease in tumor selectivity.
The model revealed that bystander effects are limited not only by diffusion and reaction of the effector, but also by effector washout, unless this washout targets downstream regions via a blood-borne bystander effect. The blood vessels act as a sink for effector, a problem that has previously been flagged in the context of drug diffusion into peritoneal tumors from the peritoneal cavity (72). Therefore, the choice of optimal effector diffusibility (defined by k tE and D E ) depends on the therapeutic context, with low diffusibility generally better for combination with radiation, and high diffusibility preferred for monotherapy activity (Figures 7B,E).
Notably, a greater increase of bystander killing with increasing effector diffusibility has previously been suggested by the finding that bystander effects of dinitrobenzamide mustard prodrugs in multicellular layer co-cultures of NfsB-overexpressing and WT cells correlated with effector lipophilicity (73,74). This may partially be due to the intimate mixture of activator and target cells in these co-cultures, so that effector washout at the multicellular layer boundaries affects both cell types to the same extent. In contrast, blood vessels in the virtual tumor microregion that act as a sink for effector are generally more distant from hypoxic prodrugactivating regions than from bystander target zones. Therefore, multicellular layer co-cultures have may limited use in assessing antitumor activity for HAP, although they are useful to rank HAP according to the diffusion range of their effectors.
Although an increase in the effector diffusion range provided only a moderate increase in single-agent activity it was found to have a surprisingly large effect in protecting cells in normal tissues where the HAP is activated (Figures 7C,F). Increasing effector diffusibility and prodrug activation rate simultaneously (the latter to compensate for washout) may therefore be a good strategy to improve monotherapy activity. The suggested HAP model parameters could be modulated in drug design by changing the physicochemical properties that determined these parameters; e.g., one-electron reduction potential in controlling k metP,max (9,75,76) or lipophilicity, pK a , molecular weight, and number of H-bond donors and acceptors in controlling diffusibility (77). In most cases any change will affect several model parameters, but the value of the SR-PK/PD model lies in part in its ability to accommodate all such changes.
The high prodrug diffusibility and low effector diffusibility that would be ideal for combination with radiotherapy according to the model predictions may not be achievable in all chemical classes of HAP. In case of dinitrobenzamide mustards (PR-104 analogs), the physicochemical properties of prodrug and effector (other than mustard reactivity) will be similar since their structures only differ in one aromatic substituent. The use of fragmenting prodrugs of aliphatic mustards such as TH-302 relaxes this constraint, providing a higher scope for independent optimization of prodrug and effector diffusion properties.

LIMITATIONS OF THE MODEL
The present SR-PK/PD model has two key limitations. Firstly, the model uses a small tumor region (0.12 mm 3 ) that cannot be expected to fully capture the heterogeneity in vascular density and Frontiers in Oncology | Pharmacology of Anti-Cancer Drugs perfusion in tumors (78). This heterogeneity results in macroregional variations in oxygenation and prodrug delivery that may decrease the overall relevance of bystander effects of Class II HAP. For example, macroscopic, well-perfused, and oxygenated areas [observed in some tumor xenografts (35,79,80)] would be beyond the reach of the hypoxia-driven bystander effects studied here, unless there is redistribution of diffusible effectors between differently oxygenated tumor regions via the bloodstream (bloodborne bystander effect). Mapped tumor regions of a larger scale are not yet available but could be used to investigate these macroregional effects in the future. Secondly, our steady-state model does not account for time-dependent processes such as cycling hypoxia (81)(82)(83) and the kinetics of reoxygenation, cell proliferation, and cell death. Cycling hypoxia is thought to arise from variations in blood flow (84)(85)(86), which would affect not only delivery of oxygen but also that of the HAP. Since cycling hypoxia has a dominant periodicity of two to three cycles per hour (82), it may be relevant even for HAP with a residence time in the circulation of few hours and critically important for HAP with long residence times in tumors as is reoxygenation for effectors with long residence times such as AQ4N (87). The implications of variable blood flow and cycling hypoxia should be incorporated in future SR-PK/PD models.

CONCLUSION
In spite of its limitations, the SR-PK/PD model has identified strategies for optimization of HAP worthy of further investigation. This is the first modeling study that attempts to systematically compare the activity of Class I and Class II HAP and the first that integrates the use of a normal tissue network to assess tumor selectivity. It suggests that class II HAP offer the following advantages relative to class I HAP: they tolerate a wider range of k metP,max (allowing a wider scope for prodrug design as well as their use in a range of human tumors with different one-electron reductase expression), are less sensitive to prodrug diffusibility (thus making it less critical to optimize this parameter) and have lower toxicity in normal tissues with physiological hypoxia (Figures 4C,F). A possible disadvantage is that complementation of radiation falls off rapidly with increasing effector diffusion away from the hypoxic activating region (Figure 7).
The modeling also showed that the design of optimized Class II HAP is complex due to the need to consider many different processes such as metabolism, membrane transfer, paracellular diffusion, and washout. The spatial and temporal heterogeneity of O 2 and blood flow within tumors (and potentially variations in reductase expression and other determinants of sensitivity) adds further complexity. Therefore it is difficult to design one HAP that fits all tumor types and therapeutic settings. However, important tendencies have been identified by the model, e.g., that antitumor activity may be increased by optimizing effector stability and prodrug activation rates while tumor selectivity may be improved by increasing effector diffusibility.
The present study has implications for targeted anticancer prodrug approaches, as most are limited by spatial heterogeneity in prodrug activation (3,4,88). The approach is also potentially applicable to hypoxia-targeting strategies that utilize prodrugs such as gene-directed enzyme prodrug therapy approaches in which the prodrug-activating enzyme is delivered using macrophages that accumulate in hypoxic tumor regions (89), or obligate anaerobic bacteria (90). In approaches with a different spatial distribution of prodrug-activating regions relative to non-activating regions and blood vessels the optimum conditions may differ and could be identified in the future using SR-PK/PD modeling.
www.frontiersin.org   5 and 7). b Ratio of overall cell killing in the tumor and cremaster muscle microregions (see Figures 5 and 7). c Overall effector released into the circulation by the FaDu tumor microregion (see Figure 8).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at http://www.frontiersin.org/Journal/10.3389/fonc.2013. 00314/abstract Figure S1 | Impact of blood flow rate on tumor PK/PD of HAP. SR-PK/PD model simulations were performed using the default O 2 transport parameters in Table 1 (Q = 40 nl/min; pO 2 = 40 mm Hg; black) and using a higher blood flow rate and lower inflow pO 2 (Q = 100 nl/min; pO 2 = 29 mm Hg; red). The prodrug activation rate constant k met P, max was set to a high value of 1 s −1 with all other parameters as in Table 2