#### Abstract

This paper presents a three-dimensional dynamic model for the chemotherapy design based on a multiphysics and multiscale approach. The model incorporates cancer cells, matrix degrading enzymes (MDEs) secreted by cancer cells, degrading extracellular matrix (ECM), and chemotherapeutic drug. Multiple mechanisms related to each component possible in chemotherapy are systematically integrated for high reliability of computational analysis of chemotherapy. Moreover, the fidelity of the estimated efficacy of chemotherapy is enhanced by atomic information associated with the diffusion characteristics of chemotherapeutic drug, which is obtained from atomic simulations. With the developed model, the invasion process of cancer cells in chemotherapy treatment is quantitatively investigated. The performed simulations suggest a substantial potential of the presented model for a reliable design technology of chemotherapy treatment.

#### 1. Introduction

Palliative chemotherapy has demonstrated improvements in survival and clinical trials for cancer patients [1]. The effectiveness of chemotherapy depends upon various factors including the type of cancer and location of cancer cells in the body [2]. The estimation of effectiveness is moreover ambiguous since chemotherapeutic drug that destroys abnormal or cancer cells also damages host cells. Because of the low efficacy of current treatments for metastatic disease, there is considerable interest in developing new systems to test chemotherapeutic drugs more quickly and accurately [3].

In the late 1980s, the National Cancer Institute (NCI) developed an in vitro drug discovery tool to test new therapeutics using 60 different human cancer cell lines [4]. The ability to test new anticancer therapeutics offers a great advantage in the development of therapies to treat cancer cells. Some reports indicated that many cancer chemotherapeutic agents are not usually selectively delivered to tumor tissues and normal cells are also destroyed by the drug [2, 5]. When administered systemically, these agents are distributed to both normal and tumor tissues by normal diffusion through capillaries [6]. An in vitro experiment has shown that the transport of small drug molecules through surrounding normal tissue is generally considered to be relatively faster than in cancer cells which are densely packed in islet [7, 8]. This finding demonstrated the reason why normal cells are destroyed faster than cancer cells by the drug. Thus, drug treatment strategy is very crucial in cancer chemotherapy and contemporary preclinical models have employed it in cancer drug development both in vitro and in vivo [5, 9]. Mathematical models of chemotherapeutic drug have been suggested to design treatment strategies that effectively destroy tumor cells while limiting toxicity on normal cells [10–14]. One of the mathematical models aimed to minimize the total amount of drug used such that the tumor cell population at the end of the treatment period reaches a specified value [11]. But it is difficult to correctly formulate and use the constraint expressing toxicity due to the cumulative drug. Recently, simulations have been performed with some proposed models [15–23]. They introduced that cancer cells had an ability to degrade and migrate into surrounding extracellular matrix (ECM). Each component involved in invasive process of tumor cells was considered as individual identity, and the model incorporated only cell-cell interaction excluding the effects of cell-matrix interaction [15]. The other mathematical models were developed to study the interaction of cancer cells with surrounding matrix [16]. But this model did not directly incorporate the interaction between cancer cells and ECM. They assumed that invasive process is triggered by cancer cells releasing matrix degrading enzymes (MDEs) which modify the ECM by degrading it. They made the effort on studying the interaction between MDEs and the degraded ECM. Byrne et al. presented a simple one-dimensional model of trophoblasts invading the uterine tissue to study placental development [17]. They explained the dominant migratory mechanism without any morphological analysis in the process of invasion. Harley et al. discussed the existence of traveling wave solutions of a haptotaxis model for malignant tumor invasion [18]. This work did not consider any energetic aspects such as surface or interface energy, which limits the accuracy of the solution. Andasari et al. studied the process of cancer-cells invasion based on mathematical analysis and computational simulation [19]. This work offered us a perspective on the ability of cancer cells to break out of tissue compartments and invade local tissue. Their simulation results were obtained by one- and two-dimensional models that limited the accurate and comprehensive understanding of the invasion process of cancer cells. In our work, a three-dimensional model incorporates multiple kinetics and energetics is developed to overcome the aforementioned limits on understanding the mechanism and morphological strategy of cancer-cells invasion. Moreover, a multiscale approach is employed, which increases the accuracy of simulation results by providing material properties with atomic simulations.

Here, we propose a three-dimensional dynamic model for the chemotherapy design. The model incorporates multiple components to study interactions between cancer cells, MDEs, degrading ECM, and chemotherapeutic drug. Multiple mechanisms possible in chemotherapy are systematically integrated for the high possibility of adequate predictions. Migration of cancer cells driven by haptotaxis, MDEs secreted by cancer cells, degradation of extracellular matrix, and diffusion of chemotherapeutic drug are analyzed simultaneously by considering multiple mechanisms in cancer-cells invasion. The velocity, number, and morphological changes of cancer cells are affected by diffusion of chemotherapeutic drug. Meanwhile, the degradation of the ECM is influenced by the diffusion of chemotherapeutic drug and MDEs which are secreted by cancer cells. Furthermore, a multiscale approach that links atomic-scale information to macroscopic properties of drugs is adopted. Molecular dynamic (MD) simulations have been carried out to examine the diffusion of chemotherapeutic drug and enhance the fidelity of the developed model. The evolving morphologies of the microstructures related to multicomponents, multiphysics, and multiscales cause computational challenges. These are addressed by employing a diffuse interface model and the reliability and effectiveness of the model are demonstrated. A series of simulations are performed to systematically investigate the effects of chemotherapeutic drugs on cancer cells and the ECM.

#### 2. Materials and Methods

##### 2.1. Multicomponents Model

A diffuse interface model has been employed to efficiently present the morphological evolution of the microstructure system, which incorporates multiple kinetics and energetics in the multicomponent system, and demonstrated its reliability and effectiveness [24–28]. We develop the diffuse interface model that incorporates a series of multicomponents, which play decisive roles in the invasion process and chemotherapy. The model incorporates multiple kinetics that are the interdiffusion of chemotherapeutic drug and cancer-cell migration induced by haptoattractant. Multiple energetics are also incorporated, which are the interface energy among the population of cancer cells, the interface energy of the ECM, and the interface energy between cancer cells and the ECM. The phase field equations with multicomponents are driven by the reduction in the total free energy of an inhomogeneous system. The multicomponents of cancer-cell invasion with chemotherapeutic drug are displayed in the schematic drawing of Figure 1. We define a field variable by the volume fraction of cancer cells, which describes the structural domain of cancer cells in the simulations. corresponds to the pure cancer cells and to other materials such as the ECM. Define a field variable by a volume fraction of the ECM. corresponds to the pure ECM and to other materials such as cancer cells, which is dimensionless. The observed morphology of the ECM from experiments will correspond to the space in the simulation results. is defined by the concentration of the MDEs and by concentration of chemotherapeutic drug. In our work, we use four variables, , , , and , which are time dependent and spatially continuous functions to describe cancer cells, MDEs, ECM, and chemotherapeutic drug. As an initial condition, we assume that corresponds to initial concentration of chemotherapeutic drug performed in the model. Various initial concentrations of chemotherapeutic drug 0.0; ; ; are tested to provide analytical information of the capability of the drug in cancer treatment. The population of cancer cells is driven by haptotaxis in the model. Haptotaxis is the directional migration of the cells due to increasing haptoattractant gradients mediated by specific material. Fibronectin is indicated to be a chemoattractant or haptoattractant for certain tumor cells [29, 30]. The maximum concentration of fibronectin is assigned on the right region of the ECM. The linear haptoattractant gradient field is introduced only along the direction to ignore the effect of the surrounding boundary on the cancer-cells invasion. Driving mechanisms for the evolution of each component as well as the interactions of them are incorporated in the following sections.

##### 2.2. Multiphysics Model

A three-dimensional dynamic model based on multiphysics technology to study the process of cancer-cells invasion with chemotherapy is expounded here. The morphological evolution of cancer cells and the ECM are represented by the free energy composed of the respective field variables, and . From the general diffuse interface approach, the total free energy of the microstructure system is followed by Cahn-Hilliard model [31] and then given byIn order to describe each component, the term could be any function with double wells. In numerical simulations, term with three components is derived as [32]where is a positive constant, a domain without cancer cells, and ECM is represented by , which is equal to . We assume for simple calculations. The remaining terms illustrate the interface energy among components. , , and are gradient energy coefficients. The net fluxes of and are given by , , and the details of the derivation of these expressions are represented in the following paragraphs.

Chemical potentials and are related to the free energy of the system and defined by . Driving forces for the evolution of each component are attributed to the chemical potential by and . We represent fluxes of and as and . Here and are mobilities of cancer cells and the ECM. The mobility of cancer cells is described by and will vanish outside the interfacial region of [24, 28, 33]. The mobility has a strong dependence on the local structure and the variable . The function of mobility instead of a constant mobility is employed and defined considering the excellent numerical stability [33], at the same time, following the model for the taxis phenomenon, which is developed by Keller and Segel [34]. The flux of is additionally induced by external stimuli and given by , where is haptotaxis sensitivity. is haptoattractant density, which represents external stimuli to provide the driving force for cancer-cells directional migration. We employed a function with a linear gradient as in our simulations. corresponds to the haptotaxis sensitivity since is proportional to the mobility of cancer cells [35], where is the sensitivity constant. The additional flux of cancer cell induced by haptotaxis can be rewritten as . Thus, the final net flux of and can be expressed as and :

The evolution of each field variable and is governed by Cahn-Hilliard nonlinear diffusion equation and combined with mass conservation relation , . The expressions of and are obtained as follows:

MDEs secreted by cancer cells and the chemotherapeutic drug are modeled as diffusing throughout the tissue and undergoing some form of decay. The kinetic process of MDEs and chemotherapeutic drug in the invasion process are described as follows:where is the mobility of MDEs and represents the released rate of MDEs by cancer cells. is the natural degradation rate of the enzymes. is the diffusion coefficient of the drug, which is obtained from the atomic calculations. The drug diffuses throughout cancer cells and the tissue and undergoes some form of decay related to .

##### 2.3. Multiscale Model

The derived governing equations are based on the continuum framework, which allows us to acquire a systematic picture of the invasive process of cancer cells and practically design the effective chemotherapy. To achieve the atomistic accuracy and broaden the applicability of the model, we employ a multiscale approach that bridges the gap between the atomic-scale information and the continuum analysis by transferring the calculated characteristics of the system with atomic simulations to the continuum calculation. Atomic simulations could provide information on calculating static and dynamic properties such as diffusion coefficient [36]. The diffusion of chemotherapeutic drug is investigated with atomic-scale calculations to enhance the possibility of expansion to various chemotherapeutic drugs.

The diffusion coefficient of chemotherapeutic drug in a polymer has been calculated from MD simulations. A cellular biological medium consists of extracellular space and extracellular polymer fibers. The amorphous polymer (polydimethylsiloxane, PDMS) is thus taken as a diffusion medium, which represents the ECM in our atomic simulations, and doxorubicin (C27H29NO11) that is chemotherapeutic drug and regarded to diffuse in PDMS. Doxorubicin is one of the most frequently used anticancer drugs in cancer treatment [37]. The diffusion coefficient is unique parameter for each specific drug and denotes the capacity of the drug diffusing in a specified medium, which is primary depending on its own chemical structure and affected by the external conditions such as the medium and temperature. Thus, the required information for calculation is the chemical structure of doxorubicin, PDMS, and temperature. Our atomic simulation is achieved with the structure of systems in a constant volume and constant density of the molecules at room temperature. Here we use a system consisting of five chemotherapeutic drug (doxorubicin) molecules and ten polydimethylsiloxane (PDMS) molecules. All the molecules are packed together and rotated. The three-dimensional drug-polymer structures are constructed by using geometry optimization method. Periodic boundary conditions are imposed in the model. Simulations are running at the temperature of 298 K with a constant volume. Then the whole structure is allowed to optimize the energy and equilibrated at the certain temperature. Purpose of the energy minimization is to ensure that dynamic simulation starts with a relatively stable structure.

Dynamic simulation begins with the optimized structure and the simulation time is selected to be 15 ps. Li et al. pointed out that value of diffusion coefficient decreased as the simulation time increased and the calculated results were consistent with the experimental value when the simulation time was between 10 ps and 20 ps [38]. In a short simulation time, the mean-square displacements show a nonlinear function with respect to the time and this tends to overestimate the diffusion coefficient [39]. With longer simulation time, the diffusion coefficient will not be changed considerably. Thus, the dynamic simulation time is selected to be 15 ps. In MD simulations, the diffusion coefficient of the drug, which is obtained from the calculated trajectory of it, can be defined as [38]where is the diffusion coefficient of the drug. is the total number of drug molecules. The angular brackets represent mean-square displacement (MSD) of the drug. Diffuse behavior can be recorded in a trajectory file for each time step of dynamic simulation. MSD obtained from the trajectory of the drug molecules, which is a time dependent function, represents an average of squared distances summed up over all possible positions of the origin. is proportional to the slope of MSD.

The diffusion coefficient of a material diffusing in a specified medium can be affected by external conditions such as concentration of the medium and temperature. We consider constant external conditions and corresponding constant diffusion parameter, which is obtained by MD simulations before the macroscopic migration simulations, to obtain the sufficient efficiency of simulations for the practical analysis. Meanwhile, thirty MD simulations with the same initial conditions have been performed to achieve the accurate average value of diffusion coefficient. Once diffusion coefficient is calculated from MD simulations, the calculated diffusion coefficient is employed to the phase field model to investigate dynamic evolution of chemotherapeutic drug. In phase field model, denotes the volume fraction of chemotherapeutic drug. As an initial condition, corresponds to initial concentration of chemotherapeutic drug considered in the model. The diffusion coefficient of chemotherapeutic drug is space and time dependent function and defined by the field variables of and , [24, 28, 33]. The diffusion coefficient of chemotherapeutic drug in the ECM is obtained from MD simulations. There is no viable result to investigate the diffusion coefficient of chemotherapeutic drug in cancer cells. Some experimental observations indicated that drug molecules diffuse with the higher diffusivity in the ECM than in cancer cells [7, 8]. Suppose that the diffusion coefficient of chemotherapeutic drug in cancer cells is approximately 50% of that in the ECM. is a dimensionless number normalized by the calculated diffusion coefficient that is obtained from MD simulations.

From the experiment, it is introduced that the treatment of solid tumors is difficult to achieve the desirable results because the chemotherapeutic drug reaching to target cells is based on diffusion [6]. The diffusion coefficient of a particular anticancer drug in a specified medium provides information about measuring whether the drug would diffuse in a favorable state. In in vitro or in vivo systems, the measurement of delivery of chemotherapeutic drug at the precise site can be facilitated by the computational results. From a targeting perspective, providing the different diffusion coefficients of chemotherapeutic drug in different tissues is useful in improving clinical treatment on tumor cells.

##### 2.4. Numerical Implementation

The process of cancer cell invasion consists of two steps; cancer cells produce MDEs and the ECM is degraded by MDEs. Meanwhile, it is well known that cancer chemotherapy is a systemic treatment. When chemotherapeutic drug is introduced, both cancer cells and the ECM are destroyed [2, 5, 6]. However, the ECM is degraded not only by the drug but also by MDEs that are secreted by cancer cells. Before the numerical discretization of the derived equations, we assume that ; here, represents the significance of the haptotaxis. Normalize the governing equations of each multicomponent with a characteristic length and time . The normalized equations of cancer cells and ECM are given bywhere and . is the degradation rate of the ECM by MDEs. The last terms of (8) and (9) denote the capacity of chemotherapeutic drug in destroying cancer cells and normal cells, respectively. The initial haptoattractant density and the mobility , , are dimensionless numbers normalized by and . The values of the parameters used in the multiscale system are listed in Table 1. The significance of the interface energy between the components is described by the Cahn numbers, . The choice of the magnitudes of the characteristic quantities depends on physical details to be resolved and computational convenience. is chosen to be 2.0 that corresponds to an experiment value of haptoattractant concentration of fibronectin, which is 200 *μ*g/ml, due to ml/*μ*g [29]. A semi-implicit Fourier spectral method is implemented here to have a high spatial resolution to resolve the high-order derivatives in the derived equations as well as the large gradients at the interface region [24]. This method is to treat the linear term implicitly and the nonlinear term explicitly to allow for larger time steps without losing numerical stability and satisfies the requirements for the numerical effectiveness [40]. Furthermore, the Semi-Implicit Backward Differentiation Formula (SBDF) scheme is applied to solve the kinetic equations without a harsh time-step constraint which has the strongest high-modal decay among the second-order multistep methods [41].

To deal with the variable mobility, the right-hand sides of (5), (8), and (9) need to be rewritten aswhere , are linear components of , . , , and . Linear terms , , and are implicitly treated and nonlinear terms , , and are explicitly treated based on the semi-implicit Fourier spectral method. There are various ways to handle the linear component [40, 41]. The numerical stability has been achieved by taking the linear term . Combined with the semi-implicit Fourier spectral method and the SBDF time integration scheme, we obtain the following discretized form: where , , , and are the constants and take , , , , andTaking the Fourier transform to (11), the governing equations of each component are given aswhere , , , , and and , , and . The caret “” and the subscript stand for the Fourier transform. , , , and can be calculated by the followed sequence. First, we compute chemical potentials that correspond to the distributions of , , , and . Then we can obtain , , , and from (13). Finally, the new , , , and are obtained from , , , and by the inverse Fourier transform. The simulation advances by repeating the procedure.

#### 3. Result and Discussion

##### 3.1. Simulation Details

Here, we quantitatively investigate the characteristics of cancer cell invasion and then systematically analyze the effect of chemotherapy on the invasion process with the developed model. The morphological evolution and migration of cancer cells and the toxicity of chemotherapeutic drug on the normal cell are specifically explored. Initial configurations of cancer cells and the ECM within the domain are illustrated in Figure 2. The domain size is set to be . The characteristic length is taken to be 3.0 *μ*m and cancer cells are initially positioned at in the spatial domain. The radius of a cancer cell is considered to be mm, which gives a typical cancer-cell volume of mm^{3} per cell [42, 43]. Initial numbers of cancer cells, and normal cells that reside in the ECM, are approximately and , respectively. Real time is equal to time scale multiplied by time steps of simulation. Time scale is chosen to be 3.2 s, due to that corresponds to an experiment value of the surface energy of poly(dimethylsiloxane) (PDMS), which is mJ/m^{2} [44]. In the following sections, we use the time step instead of real time of simulation to represent the evolution of cancer-cells invasion with chemotherapy.

##### 3.2. Diffusion of Chemotherapeutic Drug through Cancer Cells and ECM

We have performed simulations that adequately expound the diffusion process of chemotherapeutic drug through both of cancer cells and the ECM before investigating the efficacy including the systemic effect of the drug. Figure 3 shows cross-sectional views of cancer-cells invasion with chemotherapeutic drug in the computational model. The color bar presents that the field variables of cancer cells and the ECM are changing from 1.0 to 0.0. One has that (red color); (blue color); (purple color); (green color), applying with and with . The red area at the center of the domain represents cancer cells and the purple color in the surrounding area depicts the normal matrix. In the process of cancer-cells invasion, MDEs play an important role in degrading the ECM. Cancer cells secrete MDEs that destroy the normal tissue. As shown in Figure 3, MDEs are secreted by cancer cells and expressed as the varicolored region in the ECM. The number of MDEs is increasing as time increases because cancer cells secrete MDEs continuously. Chemotherapeutic drug that is described by a doughnut-shape with a dark gray color surrounds cancer cells. The mathematical model allows the drug to diffuse through both cancer cells and the ECM. Chemotherapy aims to kill abnormal or cancer cells but the normal matrix is also destroyed by the drug, which is plausibly simulated with the developed model. The clear observations supply a particularly intuitive approach and an easy way for understanding the detailed process of invasion with chemotherapy. Specialized protrusions, especially, include chemotherapeutic drug and MDEs secreted by cancer cells which are both incorporated in the dynamic model for comprehensive analysis. The three-dimensional model provides insights into understanding the mechanisms of each component in cancer-cells invasion viewing microstructural evolution and focuses our attention on specialized protrusions involved in multicomponent of cancer-cells invasion and chemotherapy treatment. Nevertheless, one- and two-dimensional models still provide important insights into cancer-cells invasion which limits our understanding of multimechanisms and different migration strategies in cancer-cells invasion. It has been possible to overcome these difficulties by creating the three-dimensional model, which presents the real state of cancer-cells invasion more precisely.

MD simulations have been carried out to study the diffusion process of chemotherapeutic drug. The diffusion coefficient can be calculated from the trajectory of the drug. MSD has been plotted as a function of time in Figure 4. The diffusion coefficient of the drug is proportional to the slope of MSD. Thirty cases have been simulated to get the accurate average value of diffusion coefficients. As shown in Figure 5, the diffusion coefficient in 60% cases is approximately cm^{2} s^{−1}. The average value of diffusion coefficient of the drug, , is cm^{2} s^{−1}. The calculated value is consistent with previous study that presented the base value of diffusion coefficient of diffusant in the cellular biological medium approximately cm^{2} s^{−1} [45]. They proposed an expression for the local effective diffusion coefficient as ; here represents diffusion coefficient in water and the base value of it is cm^{2} s^{−1}. is that is a function of the composition and fundamental geometric and physicochemical system properties, including the size of solute molecules, the size of extracellular polymer fibers, and the mass permeability of cell membrane; the base value of is mentioned and equals 0.1. They made their effort on representing the local effective diffusion coefficient of diffusant in the cellular biological medium by using the generalized function that could not provide the accurate value of diffusion coefficient of a particular diffusant in a specified medium. The other related work indicated that the diffusion coefficient of theophylline in the polymeric membrane is affected by concentration of the polymer and temperature [46]. The final diffusion coefficient of this material was cm^{2} s^{−1} at 25°C, but this value was also obtained from the numerical model. And detailed analysis of diffuse process was depending on unique structure of materials and system. Based on the above researches, it is necessary to perform detailed atomic simulation to get diffusion coefficient of the particular anticancer drug in the specified medium as accurate as possible. Our atomic simulation is primary depending on considering the detailed structure of systems and run at room temperature with constant volume and a fixed density of the polymer molecules system. We believe that the developed computational model can provide a reliable and efficient method to study the chemotherapy.

Taking the calculated diffusion coefficient of chemotherapeutic drug with MD simulation into the developed three-dimensional model of cancer-cells invasion, the diffusion process of chemotherapeutic drug is simulated as shown in Figure 6(a). The concentration of chemotherapeutic drug decreases during the chemotherapy treatment. The concentration data of the drug which is obtained from the middle region of the doughnut-shape and marked in point “” is shown in Figure 6(b) and it demonstrates that the highest concentration of the drug occurs at and the concentration is decreasing as the time increases due to natural decay phenomenon of chemotherapeutic drug. A comparison of the distribution of chemotherapeutic drug in the region “” at and is presented in Figure 6(c). It is observed that the concentration of the drug at the center region decreases and has a wider distribution as time increases since the drug diffuses into cancer cells and normal tissues. The presented simulation results adequately explain the dynamic characteristics of chemotherapeutic drug, which incorporate the atomic-scale insight obtained from atomic simulations.

**(a)**

**(b)**

**(c)**

##### 3.3. Effect of Chemotherapeutic Drug on Invasion of Cancer Cells

To quantify the effect of chemotherapeutic drug on cancer cells, we investigate the morphological changes of cancer cells without and with chemotherapeutic treatment which are shown in Figures 7(a) and 7(b), respectively. As shown in Figure 7(a), the color change from red to blue means that the variable of cancer cells is changing from 1.0 to 0.0. The change of dark red color to white color represents haptoattractant density gradient in the ECM. Higher haptoattractant density gradient in the dark red color is assigned on the right region of the ECM to induce cancer cells moving from the left to the right side. In Figure 7(b), chemotherapeutic drug is used in the system; the color change from gray to white indicates the different concentration of chemotherapeutic drug which is initially assigned as a doughnut-shape surrounding cancer cells.

**(a)**

**(b)**

As shown in Figures 7(a) and 7(b), the morphological shape of cancer cells and the number of cancer cells are changed as the chemotherapeutic drug is introduced. The population of cancer cells with chemotherapy changes its shape from spherical to deformed along the distribution of chemotherapeutic drug that hinders cancer-cells migration. As shown by two graphs, the total number of cancer cells with chemotherapeutic drug is decreased by about 14.16%. denotes the number of cancer cells. The measured velocity of cancer-cells migration is *μ*m/hr but with the presence of chemotherapeutic drug, the velocity decreases to *μ*m/hr when the time steps of simulation are . The effects of the ECM on the invasive behavior of cancer cells have been introduced by an experimental work [47]. Experimentally observed average velocity of human pancreatic cancer cells (HPAF-II) in the ECM is *μ*m/hr. The agreement is reasonably convincing.

##### 3.4. Capacity of Chemotherapy to Destroy Cells

Chemotherapy plays a role in not only affecting the morphology of the population of cancer cells but also killing cancer cells. It is also necessary to consider the effect of chemotherapy on the normal cells due to the toxicity of the drug. Many researches have made the efforts to study chemotherapy treatment in tumor cells [2, 48]. The work of delivery of chemotherapeutic agents in cancer introduced by the mathematical model revealed that approximately 60% and 80% of the initial tumor size and normal cell population reduced when the chemotherapy is performed in the treatment [2]. The seriously destroyed phenomenon on normal cells requires designing a reliable method to kill maximum number of cancer cells and minimum number of normal cells. Thus, we study the influence of chemotherapeutic drug to suggest the enhanced chemotherapy that minimizes cancer cells while limiting the toxicity on normal cells. The capacity of chemotherapeutic drug on destroying cells is represented by a parameter . The initial number of cancer cells and normal cells in the domain is considered to be and , respectively. Cancer cells killed by chemotherapeutic drug from to at are presented in Figure 8(a). As increasing the parameter , more cancer cells are destroyed by a potent drug. It is observed that at most 14.16% of cancer cells are killed by the drug when the parameter is greater than or equal to . Even though the parameter is larger than , the number of cancer cells killed by the drug does not considerably increase. Simulation results suggest the minimum efficacy of the drug with a particular amount in the cancer-cells invasion. In Figure 8(b), the reductions of normal cells with various values of parameter are presented. Approximately 11.48% of normal cells are destroyed by the drug with and the reduction of normal cells increases to 19.97% with . From the observation with Figure 8(a), the number of cancer cells is not appreciably altered with the parameter from to . The performed simulations demonstrate that the appropriate value for the capacity of chemotherapeutic drug is about in the range of to ensure that the maximum number of cancer cells and the minimum number of normal cells are destroyed under chemotherapy treatment. We consider a specific drug that corresponds to a certain diffusion coefficient and parameter . The parameter can be manipulated by a proper carrier by the development of nanotechnology [49]. Nanoparticles as a carrier provide a new mode of delivery of anticancer drug to enhance the capacity on destroying cancer cells. The motivation of the development of neoadjuvant chemotherapy is to supply an efficient and reliable drug with applicable value of the parameter on killing cells. If we consider the totally different kinds of drug that have obviously different mobility from other drugs, by recalculating the diffusion coefficient with the presented process, the expeditious analysis of chemotherapy process can be carried out with the developed model. For drugs with the recognized mobility, it would be prompt investigation. As presented by simulation results, the ECM is destroyed more seriously than cancer cells since the ECM were killed by chemotherapeutic drug and MDEs secreted by cancer cells. In addition, it is experimentally observed that drug molecules diffuse with the higher diffusivity in the ECM than in cancer cells [7, 8], because cancer cells are densely packed in islet, which leads the drug particles hard to delivery to the inside. Thus, the destroyed phenomenon on the ECM is more obvious than cancer cells, which is clearly demonstrated, and a reliable parametric study provided the developed model for the capacity of chemotherapeutic drug on killing cells.

**(a)**

**(b)**

##### 3.5. Effect of Chemotherapeutic Drug on Cancer-Cells Migration

The velocity of cancer cells as well as the destroyed amount of cancer cells could be critical for the design of chemotherapy especially when the location of cancer cells is fatal. Simulations with various initial concentrations of chemotherapeutic drug are performed to provide analytical information of the capability of the drug in cancer treatment. The velocities of cancer cells with various concentrations of chemotherapeutic drug are shown in Figure 9(a) and those at specific time steps are listed in Table 2. The maximum velocity of cancer cells is obtained with high or low concentration of drug depending on the simulation time. Before about , the maximum velocity of cancer cells is observed with the high initial concentration of drug since more ECM is degraded and more space for cancer cells migration is generated. Cancer cells, however, are degraded seriously with high initial concentration of drug. The position of cancer cells with different initial concentrations of drug is demonstrated in Figure 9(b). represents the ratio of current volume of cancer cells to the initial volume of the population of cancer cells. Within the region with the drug, cancer cells with the higher concentration of drug are more destroyed and quickly migrate through the more seriously destroyed ECM region. After when cancer cells start to escape from the region with the drug, the velocity of cancer cells with the low concentration of drug increases substantially while that with high concentration of drug decreases. Some experimental observations revealed that cancer cells could migrate through the connective tissues by their own mobility and response to the chemoattractant gradient assigned in the tissue matrix [50, 51]. The directed migration could determine that cancer cells migrate and escape from the region with the toxic effect of commonly used cancer drugs. Thus, we make the effort to incorporate this effect and investigate the velocities of cancer cells with various concentrations of chemotherapeutic drug before and after escaping from the region with the drug. Cancer cells with the high concentration of drug are degraded more seriously and secrete less MDEs than those with the low concentration of drug. Thus, after escaping from the region with the drug, cancer cells with the high concentration of drug have difficulty in generating space for invasion by degrading the ECM. Hence, the region where the drug is assigned as well as the amount of it is essential for the high efficacy of chemotherapy.

**(a)**

**(b)**

#### 4. Conclusion

In this paper, we have presented a multiphysics and multiscale model of chemotherapy and studied the invasion process of cancer cells. Cancer cells secreting MDEs, degrading ECM, and chemotherapeutic drug are systematically incorporated and the multiple mechanisms effectively integrated. Atomic simulations have been carried out and enhanced the reliability of the simulated diffusion process of chemotherapeutic drug and corresponding estimation of the efficacy of drug treatment. Simulation results demonstrate practical morphological evolution and invasive kinetics of cancer cells under chemotherapy treatment, which is limited with one- or two-dimensional models. Moreover, the simulation results performed with the developed model give insights that chemotherapeutic drug could change the morphological shape of cancer cells to prevent migration and decrease the number of cancer cells. Furthermore, the quantitative analysis provides the existence of an optimal amount of chemotherapeutic drug on improving chemotherapy treatment.

#### Conflict of Interests

The authors have no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-2013R1A1A2011263, 2014R1A2A2A09052374) and no. 2009-0083510 funded by the Korean government (MEST) through Multi-Phenomena CFD Engineering Research Center. It was also supported by Zhejiang Provincial Natural Science Foundation of China (LY15E050024) and National Natural Science Foundation of China (51175134).