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.
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 [6β8] for signaling pathway model. At the current stage in our work we have only developed the model through the blue boxes.
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).
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ββββ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 where in this system, 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 , where 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 ). In order to estimate the optimal parameters, we minimize the following objective function:
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: which is a biexponential function with four parameters through . Due to the same analysis, we adopt the same GA parameter estimation method to estimate them based on measured experimental data.
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: 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.
Also according to mass conservation law [16] we have the following: is the fluid source term from blood, and 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], was defined as 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
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 as a Neumann boundary condition for and for tumor outer boundary, where 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]: 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 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 and the outer boundary condition was set as where 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:
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:
Based on (13b) the quasi-steady-state concentration of complex can be calculated as follows:
Due to the binding loss of EGF:EGFR, the effective can be calculated by deducting drug-binded :
The effective EGF-EGFR activates down-stream effectors. High 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: where and are time-dependent tumor cell proliferation and death rate, respectively, with respect to distance to solid tumor center; is the maximum cell proliferation rate; ββis the number of occupied receptors required to generate a half-maximal response; is the maximum cell death rate; is a constant repression. An overview of the cellular signaling pathway model is illustrated in Figure 6.
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.
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 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.
(a)
(b)
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.
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).
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 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.
(a)
(b)
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.
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.
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.
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 along tumor radius at time and then compute the mean growth value , 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.
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.
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 . A longer treatment period tolerates longer treatment interval. For treatment time 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.
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).