About this Journal Submit a Manuscript Table of Contents
ISRN Biomathematics
Volume 2012 (2012), Article ID 818492, 12 pages
http://dx.doi.org/10.5402/2012/818492
Research Article

An Integrated Multiscale Mechanistic Model for Cancer Drug Therapy

1Department of Radiology, The Methodist Hospital Research Institute, Weill Cornell Medical College of Cornell University, Houston TX 77030, USA
2Department of Automation, University of Science and Technology of China, Hefei 230026, China
3Institute of Intelligent Machines, Chinese Academy of Sciences, Hefei 230031, China
4Department of Computer Science and Technology, Tongji University, Shanghai 200092, China

Received 14 August 2011; Accepted 26 September 2011

Academic Editor: A. Aouacheria

Copyright © 2012 Lei Tang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

In this paper, we established a multiscale mechanistic model for studying drug delivery, biodistribution, and therapeutic effects of cancer drug therapy in order to identify optimal treatment strategies. Due to the specific characteristics of cancer, our proposed model focuses on drug effects on malignant solid tumor and specific internal organs as well as the intratumoral and regional extracellular microenvironments. At the organ level, we quantified drug delivery based on a multicompartmental model. This model will facilitate the analysis and prediction of organ toxicity and provide important pharmacokinetic information with regard to drug clearance rates. For the analysis of intratumoral microenvironment which is directly related to blood drug concentrations and tumor properties, we constructed a drug distribution model using diffusion-convection solute transport to study temporal/spatial variations of drug concentration. With this information, our model incorporates signaling pathways for the analysis of antitumor response with drug combinations at the extracellular level. Moreover, changes in tumor size, cellular proliferation, and apoptosis induced by different drug treatment conditions are studied. Therefore, the proposed multi-scale model could be used to understand drug clinical actions, study drug therapy-antitumor effects, and potentially identify optimal combination drug therapy. Numerical simulations demonstrate the proposed system's effectiveness.

1. Introduction

Cancer disease remains a leading cause of high morbidity and mortality in both adults and children. Despite extensive efforts to discover and develop effective drugs, very few promising candidates currently exist in the development pipeline. Traditional methods for developing and testing new drugs usually involve a series of controlled experiments on groups of selected healthy volunteers to establish the relationship between dose and therapeutic effects for a given drug candidate. Traditional drug development methods are time and resource intensive and carry substantial risk for adverse effects. Hence a more efficient approach is urgently needed.

Fortunately, cancer research has undergone dramatic changes recently. One of the most important changes is the application of computational modeling for pharmacokinetics/pharmacodynamics (PK/PD) analysis [1, 2]. PK model can describe or predict the time course of drug concentration in different body compartments like blood and heart, and so forth. In contrast, PD models focus on time course of drug effects at the site of action (Figure 1). These models have been playing increasingly important roles in contemporary drug research [3]. However, if treated as two separated processes, PK and PD models may not facilitate drug research to its full extent. Some pharmaceutical companies like SimCYP and Entelos have been implementing the combined mathematical models to assist drug R&D works. The SimCYP developed a Virtual-population-based Simulator which streamlines drug development through the modelling and simulation of pharmacokinetics and pharmacodynamics analysis [4]. Entelos proposed a PhysioLab Platform for a number of immune/inflammatory dieases drug development [5]. Both of their models help reduce time/money costly drug development failures. However, none of them pays specific attention on cancer therapy and specific characters of tumor compared with normal organ. Therefore these traditional modeling approaches have limited accuracy since they did not consider the biophysical properties of different tissues and tumor interstitial space. Our approach is substantially different from their ideas. In our model, we studied how drug releases at tissue level based on the biophysical properties of tumor interstitial space. Meanwhile, we aim to propose an integrated mechanistic model for cancer drug therapy to more comprehensively evaluate drug effects at the cellular, tumor microenvironment, and organ level. To our best knowledge, this represents the first work to systematically integrate molecular imaging, biophysical modeling, and signaling pathway modeling to study how drug affects tumors in whole organisms.

818492.fig.001
Figure 1: Integrated multi-scale mechanistic model. (a) Drug effects on tumor microenvironment and different type of tumor cells via signaling pathways like EGF-EGFR, TNF𝛼-TNFR, and so forth. (b) Drug delivery from blood to tumor microenvironment.

Due to the urgency of cancer therapy and great demand to improve cancer drug development and therapeutic effects, we propose a multi-scale mechanistic model that integrates both PK and PD to study drug dose-effect relationships and to identify the optimal drug treatment strategies. The model aims at three distinct levels of detail that include the target cell, the tumor microenvironment, and whole organs where drug metabolism and/or excretion may occur. To address how drugs distribute through the whole body, we adopt a multicompartmental model to describe dynamic drug concentrations as a function of time. For unknown parameters, we refer to experimental data to make reasonable estimation which in turn may predict organ drug concentrations under different treatment conditions. With predicted blood drug concentrations, we perform a diffusion-convection solute transport model within the intratumoral space to calculate the temporal and spatial variation of drug concentration within tumor. We then use this drug distribution information as input data for the next component of the model that assesses the signaling pathway through drug binding effects. At this extracellular level, we selected the EGF-EGFR (endothelial growth factor) signaling pathway to study drug binding effects on cell proliferation and apoptosis to determine the drug treatment-therapeutic effect relationship. Our model’s detailed description of the complete mechanistic system can potentially predict the therapeutic effect and possible organ toxicity of a given drug treatment strategy, thereby enabling optimization of drug selection and dosing.

The remainder of this paper is organized as follows: Section 2 introduces the global integrated system; Section 3 describes the multi-scale model in detail, followed by Section 4, in which we present some experimental results. In the last section, some discussions and conclusions are presented.

2. Integrated System

In this section, we globally introduce the structure of an integrated system with emphasis on the function of our proposed multi-scale mechanistic model. This model aims to search for the optimal cancer therapy drug dose and treatment strategies.

The system begins at the level of macroscale drug delivery model. We use reported experimental data from molecular imaging and physiological experiments to estimate optimal parameters for the drug delivery and followed by using an optimized drug delivery system to predict concentration-time curves in different organs. Therefore, with this model we can easily change input drug dose and treatment intervals to study different therapeutic effects. With blood concentrations of drug we can further calculate spatial distribution inside a solid tumor by solving diffusion-convection partial differential equations (PDEs). Basically, drug inside tumor will accumulate and decay in accordance with drug administration and clearance. The module at the intratumoral level actually offers a way to study drug therapeutic effects by understanding how regional variations in drug distribution may affect overall antitumor effects.

The final component of the modeling approach targets intracellular mechanisms in which therapeutic response is directly linked to signaling pathway blockade or alteration. This module can provide a high-resolution overview of drug action in the targeted tumor cell. Finally, the system is designed to yield a prediction of tumor progression or positive treatment response by inputting known drug efficacy data, which includes drug dose and schedule. Hence, the first two parts in our multi-scale model which targets the organ level and intratumoral microenvironment are directly coupled to pharmacokinetics analysis which is linked to the last component of the model that analyzes the pharmacodynamic effects at the tumor cell level.

The flowchart in Figure 2 summarizes our multi-scale model at the three levels. The green dotted parallelograms indicate some potential types of data for improving this model, such as measured system coefficients for drug distribution system and reverse phase protein array (RPPA) data [68] for signaling pathway model. At the current stage in our work we have only developed the model through the blue boxes.

818492.fig.002
Figure 2: Multi-scale mechanistic model structure.

3. Model Description

The aforementioned introduce the overall structure of our proposed multi-scale mechanistic model and its primary objective for cancer therapy. Now in this part we will introduce the three modules separately to show how each modules works at different level.

3.1. Macroscale Drug Delivery Model

The first part of our model is mainly concerned with the dynamic drug concentration changes over time inside various major organs such as blood, liver, spleen, heart, and kidneys, which is essentially what constitutes PK analysis. Pharmacokinetics includes four distinct factors: drug absorption, distribution, metabolism, and excretion, which are often referred to as ADME analysis in pharmacology [9]. In our system, when studying drug absorption process, we take into consideration intravenous (IV) or oral administration since these routes represent the major forms for the vast majority of antineoplastic drugs (Figure 3).

818492.fig.003
Figure 3: Pharmacokinetics system figure. (a) Different drug administration methods including IV, oral, topical, inhalation, IM, SC, and IP. (b) Drug distribution process between blood and other organs.

For the drug distribution process, we have used a multicompartment model (Figure 3) based on a series of ordinary differential equations (ODEs) since ODEs are excellent tools for expressing dynamic system behavior. Let  [𝐶1],,[𝐶𝑖],[𝐶𝜐]  represent drug concentrations inside Organ 1 through Organ 𝜐 and let [𝐶𝐵] represent blood drug concentration; 𝑘𝑖𝐵 and𝑘𝐵𝑖 means transfer coefficients from Organ 𝑖 to blood and blood to Organ 𝑖, respectively, and so forth. Hence the drug distribution system could be depicted as𝑑𝐶1𝑑𝑡=𝑘1𝐵𝐶1+𝑘𝐵1𝐶𝐵𝑘1Loss𝐶1𝑑𝐶,(1)𝑖𝑑𝑡=𝑘𝑖𝐵𝐶𝑖+𝑘𝐵𝑖𝐶𝐵𝑘𝑖Loss𝐶𝑖,𝑑𝐶(2)𝜐𝑑𝑡=𝑘𝜐𝐵𝐶𝜐+𝑘𝐵𝜐𝐶𝐵𝑘𝜐Loss𝐶𝜐,𝑑𝐶(3)𝐵𝑘𝑑𝑡=𝐵1++𝑘𝐵𝑖++𝑘𝐵𝜐𝐶𝐵+𝑘1𝐵𝐶1++𝑘𝑖𝐵𝐶𝑖++𝑘𝜐𝐵𝐶𝜐𝑘𝐵Loss𝐶𝐵,(4) where in this system, 𝑘𝑖Loss means rate of elimination in Organ 𝑖 via metabolism or by excretion. Since elimination of a drug at equilibrium commonly follows the first-order kinetics [10], we therefore use the last term in each equation to describe drug natural decay in each organ.

In the ODE system, we need to estimate the unknown parameter, where we refer to the time course of drug biodistribution data reported in literatures. The experimental data should have a series of observations on different time points for each compartment like 𝐶𝑖(𝑡𝑗)(𝑖=1,2,𝜐,𝜐+1,𝑗=1,2,,𝑛), where (𝜐+1) is the number of compartments (including blood) and 𝑛 is the number of observations. By setting the parameters as a specific data vector Λ𝑚 (𝑚: number of parameters), we can calculate its corresponding predictions 𝐶𝑖(𝑡𝑗;Λ)(𝑖=1,2,𝜐,𝜐+1,𝑗=1,2,,𝑛). In order to estimate the optimal parameters, we minimize the following objective function:Λ=argΛmin𝜐+1𝑛𝑖=1𝑗=0𝐶𝑖(𝑡𝑗𝐶)𝑖(𝑡𝑗;Λ)2.(5)

This is a typical least square estimation objective function. But due to underlying strong noise in experimental data, some traditional least square methods (LSMs) like Gradient Descent, Gauss-Newton method, and Levenberg-Marquardt algorithm are not able to converge, yet they are prone to fall into many local optimums. We therefore use Genetic Algorithm (GA) [11] to address the global optimization problem, which is much more robust to deal with noisy data fitting problems.

When implementing this model, we firstly choose a way for drug administration, by either IV or oral. IV administration is easy to study by using the proposed ODE dynamic system since drug comes into plasma directly, either by bolus injection or by intravenous infusion; when it comes to oral administration, drug will undergo a smoother absorption process which makes it take longer time to reach peak concentration in blood and also prolongs its clearance time. We proposed the following model to fit typical oral administration blood drug concentration-time curves as shown in Figure 4:𝐶𝐵(𝑡)=𝑘1𝑒(𝑡/𝑘3)𝑘2𝑒(𝑡/𝑘4).(6) which is a biexponential function with four parameters 𝑘1 through 𝑘4. Due to the same analysis, we adopt the same GA parameter estimation method to estimate them based on measured experimental data.

818492.fig.004
Figure 4: Blood drug concentration-time curve (Oral administration).

By setting optimized parameters, the drug delivery model is able to predict drug concentration versus time curves in blood and different organs, such as 𝐶𝐵(𝑡) in blood.

3.2. Drug Distribution at Intratumoral Level

In the previous section, we discussed drug delivery at the organ level which is related to the drug concentration in blood. When drug reaches the tumor site, it will form a temporal and spatial drug distribution inside the tumor as determined by concentration changes in blood, microenvironment, and tumor properties. As solid tumors grow, development of heterogeneous neovasculature results in variable areas of rapid growth where the blood supply is adequate or necrosis at sites that are oxygen and nutrient-poor. Due to this heterogeneity, solid tumor shows heterogeneous growth patterns and has prompted particular attention on tumor angiogenesis modeling [12, 13]. Therefore our model takes into consideration of intratumoral drug distribution as a crucial role for cancer therapy study. Yet in our model we are mainly focused on the early stage of tumor treatment where the tumor is relatively small so that it could be treated as a homogeneously vascularized tiny sphere “embedded” in normal tissue. As for the physiological parameters introduced in this part, we will consider them as time independent.

3.2.1. Intratumoral Fluid Transport

As shown in Figure 5(c), intratumoral fluid convection together with diffusion effect plays an important role in delivering drug molecules. In reality, the intratumoral fluid pressure may be time dependent; however, the steady state of intratumoral fluid flow field is commonly established much faster than drug concentration field [14], especially for those macromolecule drugs. Moreover, it is assumed that the pressure field keeps unchanged without apparent changes occur of the whole tumor, like tumor radius or volume. Therefore, during the study of drug concentration, we could use Darcy’s Law to describe the quasi steady fluid transport field inside tumor [15] for simplicity:𝑢𝑖=𝐾𝑖𝜕𝑝𝑖𝜕𝑟,(7) where 𝑢𝑖 is the intratumoral fluid flow velocity, 𝐾𝑖 is the hydraulic conductivity of the intratumoral fluid, 𝑝𝑖 is the intratumoral fluid pressure, and 𝑟represents the tumor radius. Here in our model we only need to calculate the viable along the tumor radius due to spherical symmetry of early-stage tumor so as to simplify the computation.

818492.fig.005
Figure 5: (a) Cross-section of solid tumor embed in normal tissue. (b) Diffusion-convection drug transport model in solid tumor. (c) Drug transport mechanism in intratumoral space (to be continued). 𝐶𝑖 and 𝐶𝑣 are the drug concentrations in intratumoral fluid and tumor capillary respectively. Outward arrow indicates convection direction since intratumoral pressure is higher in the center than at the periphery of tumor (Figure 5(c)). (c) (Continued) Drug transport in intratumoral space.

Also according to mass conservation law [16] we have the following:𝑢𝑖=𝐹in𝐹out.(8)𝐹in is the fluid source term from blood, and 𝐹out is the drainage term mainly via the lymphatic system. Due to the lack of well-developed functional lymphatics in early-stagy tumor, we could safely neglect the drainage term in (8) as conducted in [17]. According to Starling’s hypothesis [18], 𝐹in was defined as𝐹in=𝐾𝑣𝑆𝑉𝑝𝑣𝑝𝑖𝜎𝑇𝜋𝑣𝜋𝑖.(9) Here 𝐾𝑣 is the hydraulic conductivity of tumor microvascular wall, 𝑆/𝑉 is the exchange area of the blood vessels per unit volume of the tumor tissue, 𝑝𝑣 is vascular pressure, 𝜎𝑇 is the average osmotic reflection coefficient for plasma proteins, 𝜋𝑣 is the osmotic pressure of the plasma, and 𝜋𝑖 is the osmotic pressure of intratumoral fluid. Hence we have𝐾𝑖𝑝𝑖=𝐾𝑣𝑆𝑉𝑝𝑣𝑝𝑖𝜎𝑇𝜋𝑣𝜋𝑖.(10)

There are two boundary conditions for solving this equation. For the outer boundary, we assume that the tumor tissue pressure equals the normal tissue pressure. Since the tumor has a nonflux center due to its spatial sphere symmetry, we set 𝑝𝑖|𝑟=0=0 as a Neumann boundary condition for 𝑟=0 and 𝑝𝑖|𝑟=𝑅=𝑝tissue for tumor outer boundary, where 𝑝tissue is the fluid pressure in normal tissue. We first solve (10) for intratumoral pressure distribution and then work out intratumoral fluid convection velocity according to (7). They are both important factors for drug transport at tumor site.

3.2.2. Intratumoral Drug Transport

After we have studied the intratumoral fluid, we are ready to address drug transport based on the information acquired previously. Since the heterogeneous distribution of drug is mainly induced by drug diffusion and intratumoral fluid convection, we can use convection-diffusion equation to study the drug transport process. As shown in Figure 5(c), intratumoral fluid flow transports drug molecule from a high-pressure to a low-pressure region, and diffusion effect occurs where a drug concentration gradient exits, which means that drug molecules will diffuse from high- to low-density region. For the internalization effect, we mainly refer to the chemical reaction of drug with tissue, which could be described by a first-order effective elimination process [14]. Therefore, the entire drug transport can be found in [17]:𝜕𝐶𝑖𝜕𝜕𝑡=𝐷𝜕𝑟𝜕𝐶𝑖𝜕𝑢𝜕𝑟𝑖𝐶𝑖𝜕𝑟+𝐹𝑠𝐹𝑙𝑠𝐶𝑅𝑖,(11) where 𝑢𝑖 denotes intratumoral fluid velocity as in (8), and 𝐹𝑠 and 𝐹𝑙𝑠 denote the supply and drainage term from blood and lymphatics, respectively. 𝑅(𝐶𝑖) is a sink term for drug-tissue reaction. In this equation we drop 𝐹𝑙𝑠 term as well. Also here we have𝐹𝑆=𝐹in(1𝜎)𝐶𝑣𝑃(𝑡)+𝑣𝐶𝑆/𝑉𝑣(𝑡)𝐶𝑖𝑃𝑒𝑣𝑒𝑃𝑒𝑣,𝑃1𝑒𝑣=𝐹in(1𝜎),𝑅𝐶𝑃𝑆/𝑉𝑖=𝛾𝐶𝑖,(12) where 𝑃𝑒𝑣is the ratio of convection to diffusion magnitude across tumor capillary wall; 𝑃𝑣 is the vascular permeability coefficient; 𝜎 is the osmotic reflection coefficient for the drug molecules; 𝛾 is effective reaction rate; 𝐶𝑣(𝑡) is time course of drug concentration inside tumor capillary, we assume to be the same as our predicted plasma drug concentration 𝐶𝐵(𝑡) in our macroscale drug delivery system. This time-dependent drug concentration in capillary will cause drug accumulation or decay in tumor site. For the two boundary conditions in the computational domain, due to the solid tumor’s sphere symmetry, we set the nonflux center boundary condition as 𝐶𝑖|𝑟=0=0 and the outer boundary condition was set as 𝐶𝑖|𝑟=𝑅=𝐶tissue where 𝐶tissue is drug concentration in normal tissue which could be predicted by our drug delivery system or by direct measurement. By plugging in the intratumoral fluid pressure and velocity in (7), we solve (11) for time-dependent spatial drug distribution. This drug distribution information is to be used as an input for next stage drug effect analysis.

3.2.3. Drug Effects on Tumor Cells

We will assess drug effects on a targeted tumor cell in this section. In the previous literature [19], Eladdadi and Isaacson built a mathematical model to quantitatively study the relationship between growth factor EGF, EGFR receptor numbers, and cell proliferation. An assumption was made that cell proliferation rate is a function of the number of growth factor receptors on cell surface. We take into consideration EGF-EGFR signaling pathway to represent a generic growth factor family [20], which performs as a control factor of tumor cell proliferation or apoptosis.

Normally, binding of EGF to EGFR receptor triggers complex signaling cascade with interactions between activated receptors, recruited proteins, and plasma membrane molecules, ultimately resulting in multiple down-stream effects that control cell proliferation and survival. However, in the presence of a targeted drug, this EGF-EGFR signaling pathway will be blocked and ultimately resulting in cell death in a drug dose-dependent manner. This process could be depicted as follows:EGF+EGFR𝑘𝑓𝑘𝑟EGFEGFR,(13a)EGF+EGFR+DRUG𝑘𝑓𝑘𝑟EGFEGFRDRUG.(13b)

We use the steady state of (14) to obtain the concentration of EGF : EGFR. Using Michaelis-Menten kinetics, we derive the quasi-steady state of EGF-EGFR in (13a) as follows:[]=[]EGFEGFREGFR0[]EGF[]EGF+K𝑚1.(14)

Based on (13b) the quasi-steady-state concentration of EGFEGFRDRUG complex can be calculated as follows:[]=[]EGFEGFRDRUGEGFEGFR][DRUG[]DRUG+K𝑚2.(15)

Due to the binding loss of EGF:EGFR, the effective [EGFEGFR]e can be calculated by deducting drug-binded EGFEGFRDRUG:[]EGFEGFRe=[][]EGFEGFREGFEGFRDRUG.(16)

The effective EGF-EGFR activates down-stream effectors. High [EGFEGFR]e will trigger tumor cell proliferation and reduce cell apoptosis, and vice versa. Since the number of growth factor receptors on cell surface is finite, the binding effect of drug molecule and growth factor is a nonlinear process, which means that binding rate grows fast at the beginning, yet gradually slows down when it approaches a plateau. Therefore, we depict this binding kinetics and drug therapeutic effects by using the Hill equations at Order Two:𝜔(𝑟,𝑡)=𝜔max[]EGFEGFR2e[]EGFEGFR2e+𝜇half2,(17)𝑑(𝑟,𝑡)=𝑑max1[]EGFEGFR2e/kd+1,(18) where 𝜔(𝑟,𝑡) and 𝑑(𝑟,𝑡) are time-dependent tumor cell proliferation and death rate, respectively, with respect to distance to solid tumor center; 𝜔max is the maximum cell proliferation rate; 𝜇half  is the number of occupied receptors required to generate a half-maximal response; 𝑑max is the maximum cell death rate; 𝑘𝑑 is a constant repression. An overview of the cellular signaling pathway model is illustrated in Figure 6.

818492.fig.006
Figure 6: EGF-EGFR signaling pathway and drug effects on cell proliferation and apoptosis rates.

4. Experiments

In order to show the effectiveness and advantages of our proposed model, we performed a series of experiments to illustrate the properties of our model. We analyzed both IV method and oral drug treatment since they are the most widely used administration routes. The proposed multi-scale model gives us a detailed description at different levels of drug action on the whole body and particularly in solid tumor.

4.1. Drug Delivery at Organ Level

To predict drug concentration changes in different organs, we built up the drug delivery system. For IV administration study, the system performance was to a high extent determined by system parameters introduced in (5). Therefore at first, we need to optimize the parameter set Λ𝑚 (𝑚: number of parameters) based on real-world experimental data. First we analyzed experimental data for an intravenous drug from Miao et al. [21] for an IV compound targeting on integrin pathway which belongs to the EGFR family. During the experiment, a radioactively labeled peptide, the Ecballium elaterium trypsin inhibitor (EETI-II) was synthesized and used to image tumors grown in immunodeficient mice with positron emission tomography (PET). Noninvasive time course distribution data in several organs such as the lungs, muscle, liver, tumor, and blood were acquired in their report.

In this work, we selected four compartmental data to test the model, namely, the liver, lung, tumor, and blood. After the drug delivery ODE system was set up according to (1) through (4), Genetic Algorithm was used to solve the least square objective function in (5). By setting parameters in ODE system as GA yielded, we generated the prediction data and plot them along with real data in Figure 7.

818492.fig.007
Figure 7: Experimental data fitting by Drug delivery ODE system.

Figure 8 shows drug concentration curves after bolus injection which showed rapid clearance with a half life time of approximately 2 minutes and controlled rate intravenous infusions. In clinical practice, blood drug concentrations are maintained at a constant level by intravenous infusions. With constant blood drug concentration, other organs may undergo a drug accumulation process reaching a plateau. In this figure, drug concentration rises and decays more slowly in tumor than in other organs which is consistent with the deranged neovasculature and functional lymphatics typically observed within solid tumors. Based on this knowledge, we therefore can safely eliminate the drainage terms 𝐹out and 𝐹𝑙𝑠 in (8) and (11) when solving the two equations. Moreover, as shown in Figure 8(b), the highest value of drug concentration in each organ is proportional to blood drug concentration. This relationship provides us a way to infer the maximum safe dose in blood based on maximum toxicity tolerance of different organs.

fig8
Figure 8: IV administration drug concentration curves. (a) Bolus injection. (b) Intravenous infusion.

However, if we choose drugs that are orally administration, concentration in blood will rise smoothly and fall more gradually. It takes hours to reach peak value and even longer to be cleared. This process could be described by a biexponential function model introduced in (6). For estimating the four parameters, we use the experimental data reported in [22], in which the authors measured drug concentrations in human blood when treated orally at multiple time points. Due to the possible large variations in the measured data, we also used Genetic Algorithm to estimate parameters and then plot the generated fitting curve and original data in Figure 9.

818492.fig.009
Figure 9: (Oral treatment) Experimental measured drug concentration curve and the fitting curve predicted by proposed model.

Moreover, when treated with an oral drug, patients usually need to take a dose every several hours to maintain an adequate drug concentration. The proposed oral treatment model could be used to predict drug concentration curves in plasma when treated with various timing intervals. For example, with a drug treatment interval of 12 hours, drug concentration in blood gradually rises and oscillates at a predictable level (Figure 10).

818492.fig.0010
Figure 10: Blood drug concentration curve of a fixed time interval dosing.
4.2. Drug Distribution at the Intratumoral Level

With the blood concentrations, drug distribution inside tumors can be readily calculated. This can be achieved by integrating the intratumoral fluid pressure and velocity in (10) and (7). Finite difference method (FDM) was adopted to solve the PDEs. Since tumor hosting normal tissue pressure is close to zero [23], we therefore set tumor outer boundary condition as 𝑝𝑖|𝑟=𝑅=0 for (10). Based on a diffusion-convection mechanism, pressure distribution will determine how drug diffusion from tumor capillary to tumor tissue is affected. Velocity will contribute to convection process which together with diffusion effects causes drug heterogeneous distribution. Figure 11 shows the pressure distribution (A) and velocity distribution (B) by setting equation parameters as in Table 1. We could find that the pressure is higher in the center and lower at the boundary; meanwhile velocity gradually increases. Pressure and radius are normalized in the figure.

tab1
Table 1: Baseline parameter values used for intratumoral PDE system.
fig11
Figure 11: Intratumoral fluid pressure and velocity distribution.

After obtaining pressure and velocity values, we can calculate the drug concentration based on drug transport equation (11). Assuming that tumor is growing in liver, drug concentration in liver could be used as outer boundary condition. When using intravenous infusion treatment to inhibit tumor growth, we can obtain a relationship that shows how drug inside tumor accumulates and clears (Figure 12). This approach shows that drug concentration is higher in central portions of a tumor and decreases toward tumor periphery.

818492.fig.0012
Figure 12: Intermittent intravenous infusion treatment. The treatment time is 2 hrs and treatment interval is 3 hrs.
4.3. Drug Effects at Extracellular Level

Drug distribution results provide detailed information for predicting and studying drug effects on different region inside a tumor. In this part, we use the EGFR signaling pathway model to evaluate relationship between drug concentration and treatment effects. The first part of signaling pathway model is focused on drug binding effects through the EGF receptor. Literatures clearly demonstrated that the EGF-EGFR signaling pathway activates tumor cell proliferation and specific blockade of this signaling results in the tumor inhibition [24, 25]. If drug molecules bind to EGFR so as to block EGF-EGFR signaling pathway, tumor cell proliferation will be reduced. This process was described using a series of equations from (13a) and (13b) through (18). The parameters settings are listed in Table 2.

tab2
Table 2: Baseline parameter values used for signaling pathway model.

Based on (13a) and (13b) through (18), we can calculate the drug concentration-inhibition effect relationship figure as shown in Figure 13. In the absence of drug, the ratio between tumor cell proliferation and apoptosis rate favors tumor progression. However, when drug enters the intratumoral space and binds to cell surface EGF receptor, tumor cell proliferation rate is predicted to decrease as a function of drug concentration.

818492.fig.0013
Figure 13: Drug concentration-cell growth rate relationship.

This process is described by the Hill function, and based on this model, we can calculate net cell growth rate 𝐺(𝑟,𝑐𝑟(𝑡)) at each spatial point along tumor radius by deducting cell death rate 𝐷(𝑟,𝑐𝑟(𝑡)) from cell proliferation 𝑃(𝑟,𝑐𝑟(𝑡)) rate, where 𝑟 indicates the distance from the spatial point to tumor center, and 𝑐𝑟(𝑡) is drug concentration at time𝑡. According to the analysis in the last two parts of proposed model, a specific drug dose 𝐶 will cause a particular 𝑐𝑟(𝑡) pattern.

4.4. Numerical Simulation of the Integrated Multiscale Mechanistic System

When studying different administration methods for cancer treatment, we first generated blood drug concentration curves such as that shown in Figures 8 and 10. This result was plugged into intratumoral PDE system to quantify drug distribution at the tumor site. Based on (11), heterogeneous drug gradients exist within a tumor, where drug concentration in the center is higher than that at the tumor boundary which indicates that cells at the outer boundary are more prone to divide than the other cells inside tumor. In our extracellular pharmacodynamics model, we first calculate net cell growth value 𝑅0𝐺(𝑟,𝑐𝑟(𝑡)) along tumor radius at time𝑡 and then compute the mean growth value 𝐺𝑚(𝐶)=1/𝑇𝑡𝑡+𝑇𝑅0𝐺(𝑟,𝑐(𝑡)), where 𝑇 is the length of one treatment period, and 𝐶 is the related drug dose. A typical treatment period includes treatment time and dosing interval (Figure 8) for IV or dosing interval only for oral administration (Figure 10). To inhibit tumor growing, we aimed to achieve a negative 𝐺𝑚(𝐶). With this objective in mind, we tested different treatment intervals ranging from 4 to 24 hours for oral administration. The minimum oral drug peak concentrations (P in Figure 10) with respect to forming a fundamental blood concentration pattern for inhibiting tumor growth related to each treatment interval are shown in Figure 14. When interval becomes longer, a higher drug dosage is needed.

818492.fig.0014
Figure 14: Different oral drug treatment intervals and related minimum peak concentration for inhibiting tumor growth.

When it comes to IV administration, we also focus on outer boundary drug concentration curve such as that shown in Figure 12. We have tested a series of treatment strategies with intravenous infusion time that varied from 1 to 5 hours and treatment interval times 4 to 24 hours. In order to achieve negative mean 𝐺𝑚(𝐶) in one treatment period, we searches for a minimum drug concentration for each treatment strategy. This information is presented in Table 3. It is clearly shown that when treatment interval becomes longer, the related minimum drug concentration concomitantly rises. But some of our major organs like heart, brain, and liver may reach their maximum drug tolerance levels, resulting in adverse effects. Assuming maximum organ drug tolerance related blood drug concentration is 500 nM. The concentrations that are below this value (shown in green in Table 3) should be selected. Higher concentrations (shown in red) are predicted to be toxic and should be avoided.

tab3
Table 3: Minimum drug concentration for different treatment strategy*.

The results indicate that a longer treatment period can be used with relatively longer dosing intervals. As an example, for an IV administered drug, the infusion time of 5 hours and dosing interval of less than 16 hours results in a positive cell net death rate which is predicted to inhibit tumor growth. However if the dosing interval is longer than 16 hours, the calculated cell net death rate is negative, resulting in tumor progression.

In addition, Figure 15 shows the net cell death rate evaluation for each treatment strategy when given a specific drug treatment concentration 𝐶=500nM. A longer treatment period tolerates longer treatment interval. For treatment time 𝑇=5 hours, a treatment interval less than 16 hours could ensure a positive cell net death rate that inhibits tumor growth, yet if it is longer than 16 hours, cell net death rate is negative so that tumor will remain growing.

818492.fig.0015
Figure 15: Mean cell death evaluation-IV strategy relationship. (𝑇=𝑖: treatment time is 𝑖hours; 𝐼=𝑗: interval time is 𝑗 hours).

In this part, a series of experiments have been done to demonstrate the effectiveness of the proposed multi-scale model. Experimental results indicate that the proposed model is able to serve as a powerful analysis tool to find out proper strategy for drug testing or clinical treatment.

5. Conclusion

In this work we propose a multi-scale mechanistic model for cancer drug therapy study, in which we study drug action at mechanistically distinct levels that span the macroscopic effects in major organs to the subcellular alterations in signaling cascades within tumor cells. This multiscale model provides a detailed and comprehensive view of drug action and provides important predictive relationships that may be used to cost-effectively guide both preclinic and clinical trials. However, our model is still at its early stage and several aspects must be further refined. For example, drug mechanism models at the organ level should also incorporate specific properties of drug such as composition, specific targets, and mode of action. In our initial work the ideal tumor was treated as a homogeneous tiny sphere; however tumors in practice are highly heterogeneous in size, shape, vascularity, and mode/site of invasion; these factors pose a significant challenge for accurate modeling. As already alluded, tumor growth, regional variability, and response to varying drug concentrations within tumors remain a challenge to model and to accurately determine experimentally. However such information perhaps obtainable by molecular imaging techniques may facilitate further refinements to the model. We have used a well-studied signaling pathway for the PD module in our model; however further improvement awaits greater detail in our understanding of multiple biochemical and signaling pathways involved in tumorigenesis. Nonetheless there is sufficient pharmacological and biochemical information at this time to provide a reasonably detailed framework for a comprehensive model as we have presented. Our approach provides an important foundation for analyzing and predicting drug action and antitumor effects as well as systemic toxicity to facilitate the drug discovery and development process.

Acknowledgments

The author thank Dr. Brian E. O’Neill for his kind help and suggestion and Drs. Zhengzheng Shi, Zheng Li, and Tao Peng for their discussion. This work was partially supported by the NIH Grants 5R01EB009009-03 (Li), 1R01LM010185-01 (Zhou), and NSFC Grants 61133010, 31071168, 61005010, 60905023 & 60975005 (Huang).

References

  1. P. Chaikin, G. R. Rhodes, R. Bruno, S. Rohatagi, and C. Natarajan, “Pharmacokinetics/pharmacodynamics in drug development: an industrial perspective,” Journal of Clinical Pharmacology, vol. 40, no. 12, pp. 1428–1438, 2000. View at Scopus
  2. J. Dingemanse and S. Appel-Dingemanse, “Integrated pharmacokinetics and pharmacodynamics in drug development,” Clinical Pharmacokinetics, vol. 46, no. 9, pp. 713–737, 2007.
  3. J. Y. Chien, S. Friedrich, M. A. Heathman, D. P. de Alwis, and V. Sinha, “Pharmacokinetics/pharmacodynamics and the stages of drug development: role of modeling and simulation,” AAPS Journal, vol. 7, no. 3, pp. E544–E559, 2005. View at Scopus
  4. G. L. Dickinson, S. Rezaee, N. J. Proctor, M. S. Lennard, G. T. Tucker, and A. Rostami-Hodjegan, “Incorporating in vitro information on drug metabolism into clinical trial simulations to assess the effect of CYP2D6 polymorphism on pharmacokinetics and pharmacodynamics: dextromethorphan as a model application,” Journal of Clinical Pharmacology, vol. 47, no. 2, pp. 175–186, 2007. View at Publisher · View at Google Scholar · View at Scopus
  5. A. Bangs, “Predictive biosimulation and virtual patients in pharmaceutical R and D,” Studies in Health Technology and Informatics, vol. 111, pp. 37–42, 2005. View at Scopus
  6. U. Korf, C. Lobke, O. Sahin et al., “Reverse-phase protein arrays for application-orientated cancer research,” Proteomics—Clinical Applications, vol. 3, no. 10, pp. 1140–1150, Oct 2009. View at Publisher · View at Google Scholar · View at Scopus
  7. H. Peng, J. Wen, H. Li, J. Chang, and X. Zhou, “Drug inhibition profile prediction for NFκB pathway in multiple myeloma,” PLoS ONE, vol. 6, no. 3, article e14750, 2011. View at Publisher · View at Google Scholar
  8. H. M. Peng, J. G. Wen, C. C. Chang, and X. B. Zhou, “Systematic modeling study on mechanism of p38 MAPK activation in MDS,” in The Eighteenth International Conference on Intelligent Systems for Molecular Biology (ISMB 2010), Boston,USA, 2010.
  9. A. P. Li, “Screening for human ADME/Tox drug properties in drug discovery,” Drug Discovery Today, vol. 6, no. 7, pp. 357–366, 2001.
  10. S. A. Roberts, “Drug metabolism and pharmacokinetics in drug discovery,” Current Opinion in Drug Discovery and Development, vol. 6, no. 1, pp. 66–80, 2003. View at Scopus
  11. S. Ramat and G. Magenes, “Latency detection in motor responses: a model-based approach with genetic algorithm optimization,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 10, pp. 2015–2023, 2006. View at Publisher · View at Google Scholar · View at Scopus
  12. A. R. Anderson and M. A. Chaplain, “Continuous and discrete mathematical models of tumor-induced angiogenesis,” Bulletin of Mathematical Biology, vol. 60, no. 5, pp. 857–899, 1998. View at Publisher · View at Google Scholar · View at Scopus
  13. M. A. Chaplain, S. R. McDougall, and A. R. Anderson, “Mathematical modeling of tumor-induced angiogenesis,” Annual Review of Biomedical Engineering, vol. 8, pp. 233–257, 2006. View at Publisher · View at Google Scholar
  14. C. C. Wang, J. Li, C. S. Teo, and T. Lee, “The delivery of BCNU to brain tumors,” Journal of Controlled Release, vol. 61, no. 1-2, pp. 21–41, 1999. View at Publisher · View at Google Scholar
  15. E. P. Salathe and K. N. An, “A mathematical analysis of fluid movement across capillary walls,” Microvascular Research, vol. 11, no. 1, pp. 1–23, 1976.
  16. R. K. Jain and L. T. Baxter, “Mechanisms of heterogeneous distribution of monoclonal antibodies and other macromolecules in tumors: significance of elevated interstitial pressure,” Cancer Research, vol. 48, no. 24 I, pp. 7022–7032, 1988. View at Scopus
  17. L. T. Baxter and R. K. Jain, “Transport of fluid and macromolecules in tumors. I. Role of interstitial pressure and convection,” Microvascular Research, vol. 37, no. 1, pp. 77–104, 1989.
  18. F. R. Curry, “Permeability measurements in an individually perfused capillary: the "squid axon" of the microcirculation,” Experimental Physiology, vol. 93, no. 4, pp. 444–446, 2008. View at Publisher · View at Google Scholar
  19. A. Eladdadi and D. Isaacson, “A mathematical model for the effects of HER2 overexpression on cell proliferation in breast cancer,” Bulletin of Mathematical Biology, vol. 70, no. 6, pp. 1707–1729, 2008. View at Publisher · View at Google Scholar · View at Scopus
  20. X. Zhu, X. Zhou, M. T. Lewis, L. Xia, and S. Wong, “Cancer stem cell, niche and EGFR decide tumor development and treatment response: a bio-computational simulation study,” Journal of Theoretical Biology, vol. 269, no. 1, pp. 138–149, 2011. View at Publisher · View at Google Scholar
  21. Z. Miao, G. Ren, H. Liu et al., “An engineered knottin peptide labeled with 18F for PET imaging of integrin expression,” Bioconjugate Chemistry, vol. 20, no. 12, pp. 2342–2347, 2009. View at Publisher · View at Google Scholar · View at Scopus
  22. P. A. Vasey, M. Gore, R. Wilson et al., “A phase Ib trial of docetaxel, carboplatin and erlotinib in ovarian, fallopian tube and primary peritoneal cancers,” British Journal of Cancer, vol. 98, no. 11, pp. 1774–1780, 2008. View at Publisher · View at Google Scholar · View at Scopus
  23. M. Milosevic, A. Fyles, D. Hedley et al., “Interstitial fluid pressure predicts survival in patients with cervix cancer independent of clinical prognostic factors and tumor oxygen measurements,” Cancer Research, vol. 61, no. 17, pp. 6400–6405, 2001. View at Scopus
  24. H. W. Lo, S. C. Hsu, and M. C. Hung, “EGFR signaling pathway in breast cancers: from traditional signal transduction to direct nuclear translocalization,” Breast Cancer Research and Treatment, vol. 95, no. 3, pp. 211–218, 2006. View at Publisher · View at Google Scholar
  25. S. Pennock and Z. Wang, “Stimulation of cell proliferation by endosomal epidermal growth factor receptor as revealed through two distinct phases of signaling,” Molecular and Cellular Biology, vol. 23, no. 16, pp. 5803–5815, 2003. View at Publisher · View at Google Scholar · View at Scopus