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

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 Z_c = −5. Nevertheless, the model has been tested using two additional values for cancer cell charge number, eventually accounting for different malignant activity, specifically Z_c = −3 and Z_c = −7. Glycoproteins, such as VN and PAI-1, are constituted by negatively charged oligosaccharides [25,26], for which we assume a ”conventional” charge number Z_v = −1 and Z_p = −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 i^th 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 = −Z_iu_ic_iE, where Z_i denotes the charge number, u_i expressed in m²·s⁻¹·V⁻¹ is the ionic mobility, c_i expressed in mol·m⁻³ is the concentration, and E expressed in V·m⁻¹ 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 i^th species, each ionic mobility is related to the proper diffusion parameter D_i as:

$$u_i=D_i\frac{F}{RT}$$

(1)

In Equation (1), F = 96,485.33 A·s·mol⁻¹ is the Faraday’s constant; R = 8.31 J·mol⁻¹·K⁻¹ and T are the molar gas constant and the absolute temperature measured in Kelvin, respectively.

Let Ω ⊂ ℝ³ 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 = (c₁ … c_n), where c_i= c_i (x,t), for i = 1…n is the concentration of the generic species in Ω, φ_i(x,t) is the flux of c_i through S, and s_i(c) is a source term taking into account contributions coming from all the interacting species in Ω, from the mass conservation we obtain the following:

$$\frac{d}{dt}{\displaystyle \underset{V}{\int }{c}_i(\mathbf{x},t)dx=-{\displaystyle \underset{S}{\int }{\mathbf{\phi }}_i(x,t)\cdot d\mathit{S}+{\displaystyle \underset{V}{\int }{s}_i(\mathbf{c})}d\mathit{V}}}$$

(2)

The divergence theorem applied to Equation (2) gives us the evolution equation for the i^th species as follows:

$$\frac{\partial ci}{\partial t}=-\nabla \cdot {\mathbf{\phi }}_i+s_i(\mathbf{c}),\hspace{1em}i=1,\dots n$$

(3)

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]:

$$\frac{\partial c}{\partial t}=\nabla \cdot \left[D_c\nabla c-c\left({\chi }_v\nabla v+{\chi }_u\nabla u+{\chi }_p\nabla p\right)-Z_c{u}_{c}cE\right]+{\mu }_{1}c\left(1-c/c_0\right)$$

(4)

$$\frac{\partial v}{\partial t}={\varphi }_{21}up+{\mu }_{2}v\left(1-v/v_0\right)-{\varphi }_{22}vp-\delta vm$$

(5)

$$\frac{\partial u}{\partial t}=\nabla \cdot \left(D_u{\nabla }_u\right)+{\alpha }_{31}c-{\varphi }_{31}pu-{\varphi }_{33}cu$$

(6)

$$\frac{\partial p}{\partial t}=\nabla \cdot \left(D_p{\nabla }_p-Z_p{u}_p{p}_E\right)+{\alpha }_{41}m-{\varphi }_{41}pu-{\varphi }_{42}pv$$

(7)

$$\frac{\partial m}{\partial t}=\nabla \cdot \left(D_m{\nabla }_m\right)+{\varphi }_{52}pv+{\varphi }_{53}cu-{\varphi }_{54}m$$

(8)

The constants c₀ and v₀ 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 D_c 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 (μ₁-rated). The new term Z_cu_ccE accounts for the interaction of the negatively charged cancer cells with the uniform electric field E, in which Z_c is the average charge number previously introduced as an ansatz and u_c 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 ϕ₂₁, and through a logistic term at a rate μ₂; it is also neutralized/degraded by PAI-1–plasmin interaction at a rate ϕ₂₂/δ, 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, D_u, α₃₁, ϕ₃₁, and ϕ₃₃, respectively.

In the evolution equation for PAI-1, i.e., Equation (7), the flux term takes into account the diffusion through the coefficient D_p, 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 α₄₁, ϕ₄₁, and ϕ₄₂.

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, D_m, ϕ₅₂, ϕ₅₃ and ϕ₅₄.

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|²ε⁻¹), 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 mm³ 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⁻¹ 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, c₀ = 6.7 × 10⁷ cell cm⁻³; for VN, v₀ = 1 nM; for uPA, u₀ = 1 nM; for PAI-1, p₀ = 1 nM; and for plasmin, m₀ = 0.1 nM. Diffusivities were scaled according to the reference value D = 10⁻⁶ cm² s⁻¹, and assuming as a reference distance L = 0.1 cm, we also obtained the time scaling factor τ = L²D⁻¹ = 10⁴ s. It follows that the nondimensional parameters introduced in the model are D_c = 3.5 × 10⁻⁴, D_u = 2.5 × 10⁻³, D_p = 3.5 × 10⁻³, D_m = 4.91 × 10⁻³, χ_u = 3.05 × 10⁻², χ_p = 3.75 × 10⁻², χ_v = 2.85 × 10⁻², μ₁ = 0.25, μ₂ = 0.15, α₃₁ = 0.215, α₄₁ = 0.5, δ = 8.15, φ₁₃ = 0.1, ϕ₂₁ = 0.75, ϕ₂₂ = 0.55, ϕ₃₁ = 0.75, ϕ₃₃ = 0.3, ϕ₄₁ = 0.75, ϕ₄₂ = 0.55, ϕ₅₂ = 0.11, ϕ₅₃ = 0.75, and ϕ₅₄ = 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 V₀ = −60 mV, ionic mobility obtained through Equation (1) was scaled according to V₀D⁻¹, while the applied electric field was scaled according to LV₀⁻¹. It follows that the nondimensional values for u_c, u_p, and E are, respectively, 8.1 × 10⁻⁴, 8.1 × 10⁻³, and 0.25.

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.

Description	Symbol	Unit	Dimensionless Parameter	Value
Cancer cell diffusion coefficient	Dc	cm2 s−1	DcD−1	3.5 × 10−4
uPA diffusion coefficient	Du	cm2 s−1	DuD−1	2.5 × 10−3
PAI-1 diffusion coefficient	Dp	cm2 s−1	DpD−1	3.5 × 10−3
Plasmin diffusion coefficient	Dm	cm2 s−1	DmD−1	4.91 × 10−3
uPA chemotactic coefficient	χu	cm2 s−1 nM−1	u0D−1χu	3.05 × 10−2
PAI-1 chemotactic coefficient	χp	cm2 s−1 nM−1	p0D−1χp	3.75 × 10−2
VN haptotactic coefficient	χv	cm2 s−1 nM−1	v0D−1τv	2.85 × 10−2
Cancer cell proliferation rate	μ1	s−1	τμ1	0.25
ECM proliferation rate	μ2	s−1	τμ2	0.15
uPA production rate	α	cell−1 cm3 s−1 nM	τc0u0−1α31	0.215
PAI-1 production rate	α41	s−1	τm0p0−1α41	0.5
VN degradation rate from interaction with plasmin	δ	s−1 nM−1	τm0δ	8.15
VN production rate from uPA–PAI-1 interaction	ϕ21	s−1 nM−1	τu0p0v0−1ϕ21	0.75
VN neutralization rate from interaction with PAI-1	ϕ22	s−1 nM−1	τp0ϕ22	0.55
uPA inhibition rate from interaction with PAI-1	ϕ31	s−1 nM−1	τp0ϕ31	0.75
uPA degradation rate from interaction with uPAR	ϕ33	cell−1 cm3 s−1	τc0ϕ33	0.3
PAI-1 degradation rate from interaction with uPA	ϕ41	s−1 nM−1	τu0ϕ41	0.75
PAI-1 degradation rate from interaction with VN	ϕ42	s−1 nM−1	τv0ϕ42	0.55
Plasmin production rate from PAI-1–VN interaction	ϕ52	s−1 nM−1	τp0v0m0−1ϕ52	0.11
Plasmin production rate from cancer cell–uPA interaction	ϕ53	cell−1 cm3 s−1	τu0c0m0−1ϕ53	0.75
Plasmin degradation rate	ϕ54	s−1	τϕ54	0.5
Ionic mobility for cancer cells	uc	m2 V−1 s−1	ucV0D−1	8.1 × 10−4
Ionic mobility for PAI-1	up	m2 V−1 s−1	upV0D−1	8.1 × 10−3
Electric field intensity	E	V m−1	ELV0−1	0.25

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 (x₁ = 0, y₁ = 0, z₁ = −0.5) to (x₂ = 0, y₂ = 0, z₂ = 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 Z_c = −5 representing a kind of ”mean” malignant activity, in other words, we assign, to cancer cells, a charge number Z_c = −3, or Z_c = −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., Z_c = −3 and Z_c = −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 Z_c = −3 and Z_c = −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 Z_c = −7 the invasion front reaches the upper xy boundary surface a little earlier with respect to Z_c = −5. Evidently, the simulation with Z_c = −5, shown in Figure 4b, represents the intermediate case between Z_c = −3 and Z_c = −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 Z_c 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⁻³ Vm⁻¹ along the negative z-axis direction, according to the data reported in [30], which in nondimensional form scales as E = 1.25 × 10⁻⁴. 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.