Next Article in Journal
Surrogate Modeling Approaches for Multiobjective Optimization: Methods, Taxonomy, and Results
Previous Article in Journal
A Continuous Model of Marital Relations with Stochastic Differential Equations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mathematical Model and Numerical Simulation for Electric Field Induced Cancer Cell Migration

Department of Civil, Energy, Environment and Materials Engineering, Università ‘Mediterranea’ di Reggio Calabria, Via Graziella 1, Feo di Vito, I-89122 Reggio Calabria, Italy
Math. Comput. Appl. 2021, 26(1), 4; https://doi.org/10.3390/mca26010004
Submission received: 30 November 2020 / Revised: 27 December 2020 / Accepted: 28 December 2020 / Published: 31 December 2020
(This article belongs to the Section Natural Sciences)

Abstract

:
A mathematical model describing the interaction of cancer cells with the urokinase plasminogen activation system is represented by a system of partial differential equations, in which cancer cell dynamics accounts for diffusion, chemotaxis, and haptotaxis contributions. The mutual relations between nerve fibers and tumors have been recently investigated, in particular, the role of nerves in the development of tumors, as well neurogenesis induced by cancer cells. Such mechanisms are mediated by neurotransmitters released by neurons as a consequence of electrical stimuli flowing along the nerves, and therefore electric fields can be present inside biological tissues, in particular, inside tumors. Considering cancer cells as negatively charged particles immersed in the correct biological environment and subjected to an external electric field, the effect of the latter on cancer cell dynamics is still unknown. Here, we implement a mathematical model that accounts for the interaction of cancer cells with the urokinase plasminogen activation system subjected to a uniform applied electric field, simulating the first stage of cancer cell dynamics in a three-dimensional axial symmetric domain. The obtained numerical results predict that cancer cells can be moved along a preferred direction by an applied electric field, suggesting new and interesting strategies in cancer therapy.

1. Introduction

The urokinase plasminogen activator (uPA) system is an enzymatic system that triggers proteolysis and degradation of extracellular matrix (ECM) proteins such as vitronectin (VN), which enable cancer cell proliferation and growth inside human biological tissue [1]. It is constituted by uPA and plasmin degrading enzymes, by a specific inhibitor, the plasminogen activator inhibitor type-1 (PAI-1), and by the VN [2]. Together, cancer cells and the uPA system, immersed in an aqueous and nutrient medium that constitutes the ECM, interact, while the cancer cell locomotion driven by diffusion as well chemotactic and haptotactic mechanisms takes place, the last two relying on concentration gradients of chemical species sensed by cancer cells.
After the initial cancer cell seeding, the tumor develops first in the avascular phase in which cells multiply but remain confined in a definite portion of biological tissue [3,4]. In fact, the proposed reaction–diffusion–taxis model refers to the early phase of tumor progression that begins after cancer cells somehow start to seed. This phase develops in vivo until the tumor mass reaches an average diameter of a few millimeters (about 2 mm) and is commonly known as the avascular phase as it precedes the vascular phase characterized by the tumor-induced angiogenesis [5,6,7], in which the tumor growth rate increases with the formation of new vasculature feeding the cells with nutrients [8,9]. Our model, therefore, does not provide any contribution related to tumor vasculature, because in the avascular phase nutrients are assumed to feed cancer cells through the healthy tissue vasculature surrounding the tumor. The final and most lethal phase is represented by invasion and metastasis [10], the mechanism allowing cancer cells to spread out from the original site through the vascular and lymphatic systems. The interaction of cancer cells with the uPA system triggers the degradation of the VN glycoprotein by means of the uPA serine protease degrading enzyme; plasmin is an activated serine protease catalyzing also the VN proteolysis, while PAI-1 is a glycoprotein which regulates excess of proteolysis, playing an inhibition role [2].
Recently the interaction of cancer cells with the uPA system in one and two spatial dimensions has been modeled using finite volume discretization [11], and mixed finite difference/finite element discretization [12]. Furthermore, the moving mesh partial differential equation (MMPDE) numerical technique, implemented with the finite element method (FEM), has been used to solve the partial differential equation (PDE) system arising from the mathematical modeling of the above biological system, for human tumors in the avascular phase [13,14,15]. In particular, the MMPDE numerical results, though in one spatial dimension, predict that tumor invasion develops and proliferates heterogeneously, depending on diffusion and crowding properties, as well on the presence of nutrients, while the behavior at the cancer/healthy cell interface is strongly related to malignancy [6,13].
In recent times, several studies have been carried out aimed at investigating the mutual interactions between nerves and tumors. Among others, in [16] the distribution of specific nerve fibers and their interactions in human carcinomas have been studied, finding that nerves were present inside tumors and around cells, and a regulatory role of nerve terminals in tumor growth and differentiation was also suggested. Neurotransmitters are chemicals released by nerve cells (neurons) whenever they are crossed by an appropriate electrical stimulus, and they play an important role in lymphocyte migration as tumor inhibitors, but, at the same time, they are involved in tumor progression and metastasis [17]. In Ayala et al. [18], experimental evidence of cancer-induced formation of new nerves (neurogenesis) in human prostate cancer was presented, together with the mechanism of perineural invasion, in which the cancer invasion occurred around nerves [19], which, all together, constituted a mutually beneficial system. In addition, it has recently been demonstrated that neurogenesis in mice prostate cancer induced tumor formation and metastasis [20,21]. Moreover, a mathematical model for tumor–nerve interaction was recently presented in [22], in which the results were in agreement with an experimental observation of a tumor-induced neurogenesis. Therefore, it is evident that novel nerve fibers infiltrate tumors through the neurogenesis mechanism, thus, allowing electrical stimuli to flow through the cancer cells.
Very recently, glucose metabolism in cancer cells has been experimentally investigated by [23,24], providing evidence that an abnormal rate of glycolysis, typical of cancer metabolism, induced the secretion of lactate anions which removed positive ions from the cancer cell surface, resulting in a net negative surface charge. As normal cells are neutral or slightly charged, such behavior is a hallmark of cancer cells, which can be exploited for diagnostic, and more importantly, for therapeutic purposes, with the use of charged and super-paramagnetic nanoprobes [23].
Here, we propose a three-dimensional (3D) mathematical model of the interaction of cancer cells with the uPA system, which, together, are subjected to an external uniform electric field, in which cancer cells immersed in the ECM behave as charged particles diffusing in the aqueous medium under the resulting electric force. Due to the cylindrical symmetry of the model, the simulations are performed over a two-dimensional (2D) axisymmetric domain, in order to obtain computational savings and improved numerical resolution. The model is derived from the model presented by Andasari et al. [11], adding new contributions coming from the interaction between the applied electric field and the charged particles present in the system. The early stage of cancer formation is simulated solving the resulting PDE system using the FEM, and the obtained results predict a possible way for controlling cancer cell movement using an external electric field, avoiding the use of nanoprobes or markers.

2. Materials and Methods

From an electrostatic point of view, serine proteases, such as uPA and plasmin, are charge neutral, as they contain serine (charge neutral), histidine (positive charged), and aspartic acid (negative charged) [18]. Cancer cells present a net negative charge due to glucose metabolism which can be up to 30 times with respect to normal cells [23,24]. Therefore, as an ansatz, we assigned to cancer cells a charge number Zc = −5. Nevertheless, the model has been tested using two additional values for cancer cell charge number, eventually accounting for different malignant activity, specifically Zc = −3 and Zc = −7. Glycoproteins, such as VN and PAI-1, are constituted by negatively charged oligosaccharides [25,26], for which we assume a ”conventional” charge number Zv = −1 and Zp = −1, respectively, a realistic value in light of the remarks contained in Stylianopoulos et al. [26].
We considered cancer cells, VN and PAI-1, to be charged particles in electrolytic solution. Starting from the constitutive relation for the transport of ions under electric field, being the flux φi of the ith ionic species given by the product between the ion drift velocity and the ion concentration, from foundations of electrochemistry, its final form reads as φi = −ZiuiciE, where Zi denotes the charge number, ui expressed in m2·s−1·V−1 is the ionic mobility, ci expressed in mol·m−3 is the concentration, and E expressed in V·m−1 is the applied electric field. VN is fixed in the hydrogel matrix [27,28], then, no mobility terms can be ascribed to such glycoprotein, while the ionic mobility for cancer cells and PAI-1 is calculated through the Nernst–Einstein relation, in which, referring to the ith species, each ionic mobility is related to the proper diffusion parameter Di as:
u i = D i F R T
In Equation (1), F = 96,485.33 A·s·mol−1 is the Faraday’s constant; R = 8.31 J·mol−1·K−1 and T are the molar gas constant and the absolute temperature measured in Kelvin, respectively.
Let Ω ⊂ ℝ3 be a portion of a biological tissue having bounding surface S, in which the components of the uPA system interact with a cluster of cancer cells, there initially seeded.
For t ∈ [0,T], and position x = (x,y,z) ∈ Ω, if c = (c1cn), where ci= ci (x,t), for i = 1…n is the concentration of the generic species in Ω, φi(x,t) is the flux of ci through S, and si(c) is a source term taking into account contributions coming from all the interacting species in Ω, from the mass conservation we obtain the following:
d d t V c i ( x , t ) d x = S φ i ( x , t ) d S + V s i ( c )   d V
The divergence theorem applied to Equation (2) gives us the evolution equation for the ith species as follows:
c i t = φ i + s i ( c ) , i = 1 , n
In our biological system we take into account five different species, in particular, cancer cells (c), VN (v), uPA (u), PAI-1 (p), and plasmin (m), obtaining the following PDE system to be solved for the concentration of each of the above species, while further biological details about the role played by each system component can be found in [2,11]:
c t = [ D c c c ( χ v v + χ u u + χ p p ) Z c u c c E ] + μ 1 c ( 1 c / c 0 )
v t = ϕ 21 u p + μ 2 v ( 1 v / v 0 ) ϕ 22 v p δ v m
u t = ( D u u ) + α 31 c ϕ 31 p u ϕ 33 c u
p t = ( D p p Z p u p p E ) + α 41 m ϕ 41 p u ϕ 42 p v
m t = ( D m m ) + ϕ 52 p v + ϕ 53 c u ϕ 54 m
The constants c0 and v0 in Equations (4) and (5) represent the maximum carrying capacity, respectively, for cancer cells and VN, and are detailed below. Terms in Equations (4) and (7) accounting for electric field interaction with, respectively, cancer cells and PAI-1 inhibitor, are clearly recognizable, the remaining ones in the PDE system, Equations (4)–(8), have already been detailed elsewhere [11,13,15]. Covering them quickly, in Equation (4), cancer cell movement is accounted for by cell diffusion by the Dc coefficient, also including contributions from chemotactic (χu- and χp-rated) and haptotactic (χv-rated) stimuli, from uPA and PAI-1, and VN, respectively; moreover, a logistic production term is present (μ1-rated). The new term ZcuccE accounts for the interaction of the negatively charged cancer cells with the uniform electric field E, in which Zc is the average charge number previously introduced as an ansatz and uc is the ionic mobility for cancer cells calculated through Equation (1).
VN is static [27,28], Equation (5), it is produced from the uPA–PAI-1 interaction at a rate ϕ21, and through a logistic term at a rate μ2; it is also neutralized/degraded by PAI-1–plasmin interaction at a rate ϕ22/δ, respectively. Equation (6), for uPA, accounts for diffusion, production by cancer cells, inhibition due to the interaction with PAI-1, and degradation by cancer cells at rates, Du, α31, ϕ31, and ϕ33, respectively.
In the evolution equation for PAI-1, i.e., Equation (7), the flux term takes into account the diffusion through the coefficient Dp, and also the diffusion of the negatively charged glycoprotein under the action of the electric field, as arising from Equation (1). The remaining terms are related to production of PAI-1 by plasmin, and its degradation due to the uPA–VN interaction, respectively, at rates α41, ϕ41, and ϕ42.
For plasmin, the evolution equation, i.e., Equation (8), accounts for diffusion, indirect production due to the PAI-1–VN and cancer cell–uPA interactions, and degradation, at rates, respectively, Dm, ϕ52, ϕ53 and ϕ54.
To solve the above PDE system, we impose zero flux boundary conditions, since the interacting species stay confined in the simulated domain. The initial conditions are c(x,0) = exp(−|x|2ε−1), v(x,0) = 1 − 0.5c(x,0), u(x,0) = 0.5c(x,0), p(x,0) = 0.05c(x,0), and m(x,0) = 0, where ε = 0.01 mm3 means that at t = 0 a small cancer cell cluster exists. The other conditions depend on the biological characteristics of the uPA system [11], keeping in mind that, initially, the domain is filled by VN.
We simulated the early stage of the avascular phase, under the action of an electric field with intensity E = 15 V·m−1 directed in the negative z-axis direction of a physical domain Ω = [0,r] × [0,ψ] × [−z,z] consisting of a cylinder with radius r = 0.5 mm and half-height z = 0.5 mm, where ψ = 360°. Due to the axial symmetry of the model, the domain volume is, therefore, defined by the 2D cross section in the rz-plane, R = [0,r] × [−z,z], revolved by 360 degrees in the circumferential direction, i.e., around the axis of the cylinder in the z direction. The PDE system has been solved discretizing the 2D cross section with 12,482 triangular elements and interpolating the dependent variables inside each element with quadratic shape functions, using the Galerkin’s method for the residuals of the differential equations [29], while a backward differentiation formula has been used for time stepping.
We used the model parameters taken from the literature, see [2,11] and references therein. The model was nondimensionalized introducing the following scaling quantities [15]: for cancer cells, c0 = 6.7 × 107 cell cm−3; for VN, v0 = 1 nM; for uPA, u0 = 1 nM; for PAI-1, p0 = 1 nM; and for plasmin, m0 = 0.1 nM. Diffusivities were scaled according to the reference value D = 10−6 cm2 s−1, and assuming as a reference distance L = 0.1 cm, we also obtained the time scaling factor τ = L2D−1 = 104 s. It follows that the nondimensional parameters introduced in the model are Dc = 3.5 × 10−4, Du = 2.5 × 10−3, Dp = 3.5 × 10−3, Dm = 4.91 × 10−3, χu = 3.05 × 10−2, χp = 3.75 × 10−2, χv = 2.85 × 10−2, μ1 = 0.25, μ2 = 0.15, α31 = 0.215, α41 = 0.5, δ = 8.15, φ13 = 0.1, ϕ21 = 0.75, ϕ22 = 0.55, ϕ31 = 0.75, ϕ33 = 0.3, ϕ41 = 0.75, ϕ42 = 0.55, ϕ52 = 0.11, ϕ53 = 0.75, and ϕ54 = 0.5.
To obtain a reference value for the electric quantities involved in our model, we considered that typical values for the membrane potential of cells were in the range from −50 to −70 mV [25]. Assuming, therefore, a reference electric potential V0 = −60 mV, ionic mobility obtained through Equation (1) was scaled according to V0D−1, while the applied electric field was scaled according to LV0−1. It follows that the nondimensional values for uc, up, and E are, respectively, 8.1 × 10−4, 8.1 × 10−3, and 0.25. Table 1 summarizes the model parameters used, showing only their nondimensional values together with the procedure used to obtain them, so that for each parameter its dimensional value can be inferred.

3. Results and Discussion

The computations were performed in the nondimensional 2D domain R = [0,r] × [−z,z], r∈[0,0.5], and z∈[−0.5,0.5], monitoring the dynamical evolution of the biological system with a dimensionless time step size δt = 0.1. The full axial symmetric 3D domain was recovered by a 360 degrees revolution of the obtained 2D data around the axis of the cylinder, as previously defined.
In order to obtain a benchmark for cancer cell progression in the absence of an applied electric field, we put E = 0 in Equations (4) and (7). In Figure 1, the 3D iso-surface cancer cell distribution is shown in the first column, while, in the second column, the cancer cell distribution is cut along three slices lying in the yz-planes, and then normal to the x-axis, for x = −0.38, 0, and 0.38. In all figures the cancer cell concentration is mapped in color scale. At t = 0, due to the imposed initial conditions, a cluster of cancer cells is isotropically distributed around the domain origin, as expected. At t = 8 (≈ 0.9 days), the surface delimiting the cancer/healthy cells is smooth; such behavior is unchanged up to t = 9 (≈ 1 day), when the cells growth reaches the boundary surfaces of the simulated domain, as evidenced by the central slice of the second column, while for t = 9.5 (≈ 1.1 days), because of the imposed boundary conditions, cancer cells begin to spread on the boundary surfaces.
In Figure 2, the plot of the reciprocal of the time step size versus the time step of the time-dependent solver is shown, providing the goodness of mesh convergence. The results of convergence plots relative to subsequent simulations, being quite similar, are not shown.
In Figure 3, we show the cancer evolution obtained by applying an electric field with intensity E = 0.25 in the negative z-axis direction. Until t = 8, the dynamics is not appreciably different from that shown in Figure 1. At t = 9, the invasion front has reached the upper xy boundary surface, central slice in the second column, but not yet the lower boundary surface. The side boundary surface of the domain, instead, appears just lapped by the invasion front, see the 3D iso-surface plot and the slices on the xy plane side. Increasing the time step, the anisotropy grows induced by the applied electric field along the z direction, while in the x- and y-axes direction, the tumor progression is unaffected by the electric field.
Figure 4 shows the concentration profile of cancer cells obtained along the vertical line cutting the simulated domain from (x1 = 0, y1 = 0, z1 = −0.5) to (x2 = 0, y2 = 0, z2 = 0.5), from now on referred to as the vertical cutline, i.e., along the line splitting each central slice, shown in the second columns of Figure 1 and Figure 3, into two symmetric domains with respect the z-axis. In the absence of an applied electric field, Figure 4a, at each time step, the concentration profile is symmetric with respect to z = 0. At t = 0, continuous blue line, the cells are distributed according to the initial conditions imposed. During the time evolution, the invasion front propagates symmetrically; at t = 9 (yellow line/cross) the invasion front just licks both xy boundary surfaces, which are definitively touched at t = 9.5 (purple line/diamond). Switching the electric field on, see Figure 4b, the profiles for t > 0 show a growing asymmetry according to the time evolution. In fact, already at t = 8 (red line/circle) the cancer cell profile appears to be asymmetric in both intensity and position of concentration peaks, the latter being slightly shifted towards increasing z values, i.e., towards the upper xy boundary surface located at z = 0.5 as compared with the case with E = 0. Such behavior is enhanced at the higher time steps, in particular, for t = 9, the cancer cell concentration increases close to z = 0.5, proving that the invasion front has almost reached the upper xy boundary surface, contrary to the lower one at z = −0.5. Finally, for t = 9.5, the maximum of the invasion front is on the upper xy boundary surface at z = 0.5, while it is not yet on the lower boundary surface.
As the malignant activity of cancer cells is related to their negative net charge, as a result of abnormal glucose metabolism [23,24], we hypothesize that the cancer cell charge number can assume values slightly below, or slightly above, with respect to the ansatz Zc = −5 representing a kind of ”mean” malignant activity, in other words, we assign, to cancer cells, a charge number Zc = −3, or Zc = −7, depending on whether the malignant activity is lower or higher than the average. In this respect, we have computed the dynamical evolution of our system subjected to the same electric field imposing two more charge numbers for cancer cells, i.e., Zc = −3 and Zc = −7. The behavior of the 3D iso-surface cancer cell distribution, and of their 2D slices, are similar, from a qualitative point of view, to the images shown in Figure 3, therefore, we do not show them; instead, we show the cancer cell concentration profiles along the previously defined vertical cutline, in Figure 5a,b, for Zc = −3 and Zc = −7, respectively. As it can be seen, during the time evolution, the asymmetries already observed in Figure 4b are still present, weaker in Figure 5a, and more pronounced in Figure 5b. By increasing the charge number for cancer cells, the tumor invasion front moves rigidly towards the upper confining surface along the electric field direction, and the intensity of the asymmetries also grows accordingly. In particular, comparing the profiles of Figure 4b and Figure 5b, it can be seen that with Zc = −7 the invasion front reaches the upper xy boundary surface a little earlier with respect to Zc = −5. Evidently, the simulation with Zc = −5, shown in Figure 4b, represents the intermediate case between Zc = −3 and Zc = −7.
To better inspect such asymmetries, in Figure 6 we show a one-to-one comparison between the profiles obtained with, and without, the applied electric field. Ongoing, from t = 8, Figure 6a, to t = 9.5, Figure 6c, in the absence of an electric field (continuous blue line), the profiles are symmetrical with respect to z = 0. Applying the electric field, for the Zc values −3 (red/circle), −5 (yellow/cross), and −7 (purple/diamond), the concentration profiles are gradually pushed towards higher z values, also showing an intensity asymmetry of the concentration peaks growing with time. In particular, from Figure 6, it can be deduced that, at each time step, the distance between the maxima of the concentration profiles is unchanged in both the E = 0 and E = 0.25 cases.
Since, along the nerve fibers infiltrating tumors, electrical stimuli are transmitted by the neurons, in turn triggering the mutual interactions between nerves and tumors [16,17,18,19,20,21,22], it could be hypothesized that the internal electric fields generated by the electrical stimuli flowing through the nerve fibers could also influence the locomotion of cancer cells. Nevertheless, in a very recent experimental study in vivo carried out on mice inoculated with mammary cancer cells [30], electrogram recordings within the tumors detected a bioelectric activity consisting of pulsed waves propagating with amplitude in the microvolt range. Such magnitudes are well below typical values for the membrane potential of cells, and even more so with respect to the values that could arise from the external electric field simulated here. In order to verify possible effects on the cancer cell locomotion induced by the internal electric field generated by the electric pulses propagating along the nerve fibers, we performed a simulation introducing in Equations (4) and (7) an electric field with amplitude E = 7.5 × 10−3 Vm−1 along the negative z-axis direction, according to the data reported in [30], which in nondimensional form scales as E = 1.25 × 10−4. Furthermore, to enhance the effect of the electric field on cancer cells, we chose a charge number Zc = −7. As expected, the numerical results, not shown, do not differ appreciably from those obtained in the absence of an applied electric field shown in Figure 1, thus, evidencing that the internal electric field possibly generated in peripheral tumors is irrelevant for electric-field driven cancer cell locomotion.
Such findings indicate that tumor progression subjected to an applied electric field is biased towards the field direction, more precisely, pushed by the electric field, negatively charged cancer cells tend to accumulate along the field direction and close to the appropriate domain boundary; in addition, the invasion front is rigidly moved along the field direction, but the invasion rate is nearly unchanged. Interestingly, the model proves to be sensitive to variations of the charge number of cancer cells, which is directly linked to their malignancy degree.

4. Conclusions

The numerical evidence, arising from the proposed mathematical model for cancer cell dynamics in the early stage of the avascular phase, suggests that a uniform electric field applied along a prescribed direction of the modeled biological domain could influence the tumor dynamics. The performed computations show that cancer cells grow during the time evolution with a spatially symmetric invasion. When subjected to an external uniform electric field, the invasion front is biased along the field direction, without an increase in invasion speed, but two distinct effects are observed, i.e., cancer cells tend to concentrate towards the field direction, and the invasion front translates along the same direction. Inside the tumor there can be internal electric fields generated by the electrical pulses flowing on the nerve fibers that infiltrate it, but our simulations indicate that they are not suitable for inducing a displacement of cancer cells. At present, experimental studies carried out on the same system, or computations simulating similar conditions are lacking, therefore, we cannot make any comparisons, while the present results must be validated by means of appropriate experiments. Nevertheless, they allow us to hypothesize that an external uniform electric field applied to a cancer cell progression could be capable of inducing and driving their locomotion, opening up an important issue related to therapeutic and surgical strategies effective for cancer cell targeting.

Funding

The APC was funded by the Italian Ministry of University and Research (MUR).

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Andreasen, P.A.; Egelund, R.; Petersen, H.H. The plasminogen activation system in tumor growth, invasion, and metastasis. Cell Mol. Life Sci. 2000, 57, 25–40. [Google Scholar] [CrossRef] [PubMed]
  2. Chaplain, M.A.J.; Lolas, G. Mathematical modelling of cancer cell invasion of tissue: The role of the urokinase plasminogen activation system. Math. Model Methods Appl. Sci. 2005, 15, 1685–1734. [Google Scholar] [CrossRef] [Green Version]
  3. Preziosi, L. Modelling tumour growth and progression. In Progress in Industrial Mathematics at ECMI 2002; Springer: Berlin/Heidelberg, Germany, 2004; pp. 53–66. [Google Scholar]
  4. Byrne, H.M. Modelling avascular tumor growth. In Cancer Modelling and Simulation; Chapman & Hall/CRC Press: Boca Raton, FL, USA, 2003; pp. 75–120. [Google Scholar]
  5. Chaplain, M.A.J.; McDougall, S.R.; Anderson, A.R.A. Mathematical Modeling of Tumor-Induced Angiogenesis. Annu. Rev. Biomed. Eng. 2006, 8, 233–257. [Google Scholar] [CrossRef] [Green Version]
  6. Araujo, R.P.; McElwain, D.L.S. A history of the study of solid tumour growth: The contribution of mathematical modeling. Bull. Math. Biol. 2004, 66, 1039–1091. [Google Scholar] [CrossRef] [PubMed]
  7. Folkman, J. Tumor angiogenesis. Adv. Cancer. Res. 1974, 19, 331–358. [Google Scholar] [PubMed]
  8. Folkman, J. The vascularization of tumors. Sci. Am. 1976, 234, 58–73. [Google Scholar] [CrossRef] [PubMed]
  9. Folkman, J.; Klagsbrun, M. Angiogenic factors. Science 1987, 235, 442–447. [Google Scholar] [CrossRef]
  10. Hanahan, D.; Weinberg, R.A. Hallmarks of cancer: The next generation. Cell 2011, 144, 646–674. [Google Scholar] [CrossRef] [Green Version]
  11. Andasari, V.; Gerisch, A.; Lolas, G.; South, A.P.; Chaplain, M.A.J. Mathematical modeling of cancer cell invasion of tissue: Biological insight from mathematical analysis and computational simulation. J. Math. Biol. 2011, 63, 141–171. [Google Scholar] [CrossRef]
  12. Peng, L.; Trucu, D.; Lin, P.; Thompson, A.; Chaplain, M.A.J. A multiscale mathematical model of tumour invasive growth. Bull. Math. Biol. 2017, 79, 389–429. [Google Scholar] [CrossRef] [Green Version]
  13. Amoddeo, A. Adaptive grid modelling for cancer cells in the early stage of invasion. Comput. Math. Appl. 2015, 69, 610–619. [Google Scholar] [CrossRef]
  14. Amoddeo, A. Oxygen induced effects on avascular tumour growth: A preliminary simulation using an adaptive grid algorithm. J. Phys. Conf. Ser. 2015, 633, 012088. [Google Scholar] [CrossRef]
  15. Amoddeo, A. A moving mesh study for diffusion induced effects in avascular tumour growth. Comput. Math. Appl. 2018, 75, 2508–2519. [Google Scholar] [CrossRef]
  16. Lu, S.; Zhou, Y.; Que, H.; Liu, S. Peptidergic innervation of human esophageal and cardiac carcinoma. World J. Gastroenterol. 2003, 9, 399–403. [Google Scholar] [CrossRef]
  17. Lang, K.; Bastian, P. Neurotransmitter effects on tumor cells and leukocytes. In Neuronal Activity in Tumor Tissue; Zänker, K.S., Entschladen, F., Eds.; Prog Exp Tumor Res.; Karger: Basel, Switzerland, 2007; Volume 39, pp. 99–121. [Google Scholar]
  18. Ayala, G.E.; Dai, H.; Powell, M.; Li, R.; Ding, Y.; Wheeler, T.M.; Shine, D.; Kadmon, D.; Thompson, T.; Miles, B.J.; et al. Cancer-related axonogenesis and neurogenesis in prostate cancer. Clin. Cancer Res. 2008, 14, 7593–7603. [Google Scholar] [CrossRef] [Green Version]
  19. Villers, A.; McNeal, J.; Redwine, E.; Freiha, F.; Stamey, T. The role of perineural space invasion in the local spread of prostatic adenocarcinoma. J. Urol. 1989, 142, 763–768. [Google Scholar] [CrossRef]
  20. Magnon, C.; Hall, S.J.; Lin, J.; Xue, X.; Gerber, L.; Freedland, S.J.; Frenette, P.S. Autonomic nerve development contributes to prostate cancer progression. Science 2013, 341, 1236361. [Google Scholar] [CrossRef] [Green Version]
  21. Magnon, C. Role of the autonomic nervous system in tumorigenesis and metastasis. Mol. Cell. Oncol. 2015, 2, e975643. [Google Scholar] [CrossRef] [Green Version]
  22. Lolas, G.; Bianchi, A.; Syrigos, K.N. Tumour-induced neoneurogenesis and perineural tumour growth: A mathematical approach. Sci. Rep. 2016, 6, 20684. [Google Scholar] [CrossRef]
  23. Le, W.; Chen, B.; Cui, Z.; Liu, Z.; Shi, D. Detection of cancer cells based on glycolytic-regulated surface electrical charges. Biophys. Rep. 2019, 5, 10–18. [Google Scholar] [CrossRef] [Green Version]
  24. Chen, B.; Le, W.; Wang, Y.; Li, Z.; Wang, D.; Ren, L.; Lin, L.; Cui, S.; Hu, J.J.; Hu, Y.; et al. Targeting negative surface charges of cancer cells by multifunctional nanoprobes. Theranostics 2016, 6, 1887–1898. [Google Scholar] [CrossRef]
  25. Nelson, D.L.; Cox, M.M. Lehninger Principles of Biochemistry; W.H. Freeman and Company: New York, NY, USA, 2008. [Google Scholar]
  26. Stylianopoulos, T.; Poh, M.Z.; Insin, N.; Bawendi, M.G.; Fukumura, D.; Munn, L.L.; Jain, R.K. Diffusion of particles in the extracellular matrix: The effect of repulsive electrostatic interactions. Biophys. J. 2010, 99, 1342–1349. [Google Scholar] [CrossRef] [Green Version]
  27. Lieleg, O.; Baumgartel, R.M.; Bausch, A.R. Selective filtering of particles by the extracellular matrix: An electrostatic bandpass. Biophys. J. 2009, 97, 1569–1577. [Google Scholar] [CrossRef] [Green Version]
  28. Anderson, A.R.A.; Chaplain, M.A.J. Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Biol. 1998, 60, 857–900. [Google Scholar] [CrossRef]
  29. Zienkiewicz, O.C.; Taylor, R.L. The Finite Element Method; Butterworth–Heinemann: Oxford, UK, 2002. [Google Scholar]
  30. McCallum, G.A.; Shiralkar, J.; Suciu, D.; Covarrubias, G.; Yu, J.S.; Karathanasis, E.; Durand, D.M. Chronic neural activity recorded within breast tumors. Sci. Rep. 2020, 10, 14824. [Google Scholar] [CrossRef]
Figure 1. Cancer cell distribution inside the simulated domain during the time evolution, without applied electric field. (first column) three-dimensional (3D) iso-surface plots; (second column) Slices cutting the plots of the first column along planes normal to the x direction. Snapshots taken at t = 0, t = 8, t = 9, and t = 9.5.
Figure 1. Cancer cell distribution inside the simulated domain during the time evolution, without applied electric field. (first column) three-dimensional (3D) iso-surface plots; (second column) Slices cutting the plots of the first column along planes normal to the x direction. Snapshots taken at t = 0, t = 8, t = 9, and t = 9.5.
Mca 26 00004 g001aMca 26 00004 g001b
Figure 2. Reciprocal of the time step size versus the time step. Increasing time step sizes indicate good mesh convergence.
Figure 2. Reciprocal of the time step size versus the time step. Increasing time step sizes indicate good mesh convergence.
Mca 26 00004 g002
Figure 3. Cancer cell distribution inside the simulated domain during the time evolution, applying a uniform electric field in the negative z-axis direction with intensity E = 0.25. (first column) 3D iso-surface plots; (second column) Slices cutting the plots of the first column along planes normal to the x direction. Snapshots taken at t = 0, t = 8, t = 9, and t = 9.5.
Figure 3. Cancer cell distribution inside the simulated domain during the time evolution, applying a uniform electric field in the negative z-axis direction with intensity E = 0.25. (first column) 3D iso-surface plots; (second column) Slices cutting the plots of the first column along planes normal to the x direction. Snapshots taken at t = 0, t = 8, t = 9, and t = 9.5.
Mca 26 00004 g003
Figure 4. Profiles of the cancer cell concentration along the vertical cutline (see text). (a) In the absence of an applied electric field; (b) Applying a uniform electric field in the negative z-axis direction, with intensity E = 0.25, putting Zc = −5. Snapshots are taken at t = 0 (blue), t = 8 (red/circle), t = 9 (yellow/cross), and t = 9.5 (purple/diamond).
Figure 4. Profiles of the cancer cell concentration along the vertical cutline (see text). (a) In the absence of an applied electric field; (b) Applying a uniform electric field in the negative z-axis direction, with intensity E = 0.25, putting Zc = −5. Snapshots are taken at t = 0 (blue), t = 8 (red/circle), t = 9 (yellow/cross), and t = 9.5 (purple/diamond).
Mca 26 00004 g004
Figure 5. Profiles of the cancer cell concentration along the vertical cutline, applying a uniform electric field in the negative z-axis direction, with intensity E = 0.25. (a) Zc = −3; (b) Zc = −7. Snapshots are taken at t = 0 (blue), t = 8 (red/circle), t = 9 (yellow/cross), and t = 9.5 (purple/diamond).
Figure 5. Profiles of the cancer cell concentration along the vertical cutline, applying a uniform electric field in the negative z-axis direction, with intensity E = 0.25. (a) Zc = −3; (b) Zc = −7. Snapshots are taken at t = 0 (blue), t = 8 (red/circle), t = 9 (yellow/cross), and t = 9.5 (purple/diamond).
Mca 26 00004 g005aMca 26 00004 g005b
Figure 6. A comparison of the profiles of the cancer cell concentration along the vertical cutline. (a) t = 8; (b) t = 9; and (c) t = 9.5. In each panel, the blue line refers to E = 0, the red line/circle refers to E = 0.25, Zc = −3, the yellow line/cross refers to E = 0.25, Zc = −5, and the purple line/diamond refers to E = 0.25, Zc = −7.
Figure 6. A comparison of the profiles of the cancer cell concentration along the vertical cutline. (a) t = 8; (b) t = 9; and (c) t = 9.5. In each panel, the blue line refers to E = 0, the red line/circle refers to E = 0.25, Zc = −3, the yellow line/cross refers to E = 0.25, Zc = −5, and the purple line/diamond refers to E = 0.25, Zc = −7.
Mca 26 00004 g006aMca 26 00004 g006b
Table 1. Model parameter values used in the simulations. The fourth column shows the procedure used to nondimensionalize the parameters. For reference values, please see text.
Table 1. Model parameter values used in the simulations. The fourth column shows the procedure used to nondimensionalize the parameters. For reference values, please see text.
DescriptionSymbolUnitDimensionless ParameterValue
Cancer cell diffusion coefficientDccm2 s−1DcD−13.5 × 10−4
uPA diffusion coefficientDucm2 s−1DuD−12.5 × 10−3
PAI-1 diffusion coefficientDpcm2 s−1DpD−13.5 × 10−3
Plasmin diffusion coefficientDmcm2 s−1DmD−14.91 × 10−3
uPA chemotactic coefficientχucm2 s−1 nM−1u0D−1χu3.05 × 10−2
PAI-1 chemotactic coefficientχpcm2 s−1 nM−1p0D−1χp3.75 × 10−2
VN haptotactic coefficientχvcm2 s−1 nM−1v0D−1τv2.85 × 10−2
Cancer cell proliferation rateμ1s−1τμ10.25
ECM proliferation rateμ2s−1τμ20.15
uPA production rateαcell−1 cm3 s−1 nMτc0u0−1α310.215
PAI-1 production rateα41s−1τm0p0−1α410.5
VN degradation rate from interaction with plasminδs−1 nM−1τm0δ8.15
VN production rate from uPA–PAI-1 interactionϕ21s−1 nM−1τu0p0v0−1ϕ210.75
VN neutralization rate from interaction with PAI-1ϕ22s−1 nM−1τp0ϕ220.55
uPA inhibition rate from interaction with PAI-1ϕ31s−1 nM−1τp0ϕ310.75
uPA degradation rate from interaction with uPARϕ33cell−1 cm3 s−1τc0ϕ330.3
PAI-1 degradation rate from interaction with uPAϕ41s−1 nM−1τu0ϕ410.75
PAI-1 degradation rate from interaction with VNϕ42s−1 nM−1τv0ϕ420.55
Plasmin production rate from PAI-1–VN interactionϕ52s−1 nM−1τp0v0m0−1ϕ520.11
Plasmin production rate from cancer cell–uPA interactionϕ53cell−1 cm3 s−1τu0c0m0−1ϕ530.75
Plasmin degradation rateϕ54s−1τϕ540.5
Ionic mobility for cancer cellsucm2 V−1 s−1ucV0D−18.1 × 10−4
Ionic mobility for PAI-1upm2 V−1 s−1upV0D−18.1 × 10−3
Electric field intensityEV m−1ELV0−10.25
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Amoddeo, A. Mathematical Model and Numerical Simulation for Electric Field Induced Cancer Cell Migration. Math. Comput. Appl. 2021, 26, 4. https://doi.org/10.3390/mca26010004

AMA Style

Amoddeo A. Mathematical Model and Numerical Simulation for Electric Field Induced Cancer Cell Migration. Mathematical and Computational Applications. 2021; 26(1):4. https://doi.org/10.3390/mca26010004

Chicago/Turabian Style

Amoddeo, Antonino. 2021. "Mathematical Model and Numerical Simulation for Electric Field Induced Cancer Cell Migration" Mathematical and Computational Applications 26, no. 1: 4. https://doi.org/10.3390/mca26010004

Article Metrics

Back to TopTop