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

: 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 ﬁbers 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 ﬂowing along the nerves, and therefore electric ﬁelds 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 ﬁeld, 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 ﬁeld, simulating the ﬁrst 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 ﬁeld, suggesting new and interesting strategies in cancer therapy.


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.

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 i u i c i E, where Z i denotes the charge number, u i expressed in m 2 ·s −1 ·V −1 is the ionic mobility, c i 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 i th species, each ionic mobility is related to the proper diffusion parameter D i as: 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 Ω ⊂ R 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 = (c 1 . . . 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: The divergence theorem applied to Equation (2) gives us the evolution equation for the i th species as follows: 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]: The constants c 0 and v 0 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 (µ 1 -rated). The new term Z c u c cE 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 φ 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, D u , α 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 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 α 41 , φ 41 , and φ 42 .
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 halfheight 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 0 = 6.7 × 10 7 cell cm −3 ; for VN, v 0 = 1 nM; for uPA, u 0 = 1 nM; for PAI-1, p 0 = 1 nM; and for plasmin, m 0 = 0.1 nM. Diffusivities were scaled according to the reference value D = 10 −6 cm 2 s −1 , and assuming as a reference distance L = 0.1 cm, we also 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 0 = −60 mV, ionic mobility obtained through Equation (1) was scaled according to V 0 D −1 , while the applied electric field was scaled according to LV 0 −1 . It follows that the nondimensional values for u c , u p , 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.
VN production rate from uPA-PAI-1 interaction Ionic mobility for PAI-1

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.
tem 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.
Math. Comput. Appl. 2021, 26, x FOR PEER REVIEW 7 of 14 t = 9.5 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 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. t = 9.5 (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.
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. 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 xand 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 1 = 0, y 1 = 0, z 1 = −0.5) to (x 2 = 0, y 2 = 0, z 2 = 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 Figures 1 and 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.   Figure 4 shows the concentration profile of cancer cells obtained along the vertica 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 centra slice, shown in the second columns of Figures 1 and 3, into two symmetric domains with 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 Figures 4b and 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. lution 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 Figures 4b and 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   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. 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 −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. 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.
(a) 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 .

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.