\begin{document}

\title*{GPU-based parallel simulations of the Gatenby-Gawlinski model with anisotropic, heterogeneous acid diffusion}
\titlerunning{Simulation on GPUs of anisotropic, heterogeneous acid-mediated tumor invasion}


\author{Corrado Mascia, Donato Pera and Chiara Simeoni}

\institute{Corrado Mascia \at Dipartimento di Matematica G. Castelnuovo, Sapienza Universit\`a di Roma,\\
piazzale Aldo Moro 2 - 00185 Roma - Italy\\
\email{corrado.mascia@uniroma1.it}
\and Donato Pera \at Dipartimento di Ingegneria e Scienze dell'Informazione e Matematica,\\
Universit\`a degli studi dell'Aquila, via Vetoio (snc) - 67100 Coppito (AQ) - Italy\\
\email{donato.pera@univaq.it}
\and Chiara Simeoni (corresponding author) 
\at Laboratoire de Math\'ematiques J.A. Dieudonn\'e CNRS UMR 7351,\\
Universit\'e C\^ote D'Azur, Parc Valrose - 06108 Nice Cedex 2 - France\\
\email{chiara.simeoni@univ-cotedazur.fr}}

\maketitle


\abstract{We introduce a variant of the Gatenby-Gawlinski model for acid-mediated tumor invasion, accounting for anisotropic and heterogeneous diffusion of the lactic acid across the surrounding healthy tissues. Numerical simulations are performed for two-dimensional data by employing finite volume schemes on staggered Cartesian grids, and parallel implementation through the modern CUDA GPUs technology is considered. The effectiveness of such approach is proven by reproducing biologically relevant results like the formation of propagating fronts and the emergence of an interstitial gap between normal and cancerous cells, which is driven by the pH lowering strategy and depends significantly on the diffusion rates. By means of a performance analysis of the serial and parallel execution protocols, we infer that exploiting highly parallel GPU-based computing devices allows to rehabilitate finite volume schemes on regularly-shaped meshes, together with explicit time discretization, for complex applications to interface diffusion problems of invasive processes.}
\medskip

\textbf{Keywords:} acid-mediated tumor invasion, anisotropic and heterogenous diffusion, propagating fronts, finite volume schemes, GPU-based parallel computing.

\medskip

\textbf{MSC(2010):} 65M08, 65Y05, 68W10, 35K57, 35Q92, 74E05, 74E10, 92C15.


\section{Biological context and mathematical modeling}
\label{sec:introduction}


The phenomenological framework concerning the \textit{Warburg effect} relies on the experimental results achieved by Otto Warburg~\cite{Warb1956} in the 1920s, which essentially prove that cancer cells metabolism leans on anaerobic glycolysis for adenosine triphosphate (ATP) production, regardless of the available oxygen amount, hence causing lactic acid fermentation. 
Such a behavior leads to the \textit{acid-mediated invasion hypothesis}~\cite{SmGa2008}, whose key point consists in assuming that acidification induced by an excess of lactic acid sets up a toxic microenvironment for normal cells and, consequently, favors malignant cells spreading.


From the point of view of mathematical modeling, these qualitative statements are properly framed by a system of reaction-diffusion equations known as the \textit{Gatenby-Gawlinski model}~\cite{GaGa1996}, which describes tumor cells proliferation and progression inside the local healthy tissue through the destructive effect of superfluous acidity, thus focusing on a stage of biological development in which the carcinogenesis has already occurred.

The crucial feature revealed also by analytical investigations is the emergence of propagating fronts connecting equilibrium states, but the available literature is currently restricted to almost only one-dimensional problems~\cite{FaHe2009,McGa2014,DaHe2018}.


Multidimensional configurations have been recently considered in~\cite{MoSi2019} to support experimental observations, by performing two- and three-dimensional finite element simulations exploiting the COMSOL Multiphysics$^{\textregistered}$ environment~\cite{CoMu2020} together with ParaView utilities~\cite{PaVi2020} for post-processing three-dimensional graphical results.

In this article, we introduce a biologically relevant variant of the Gatenby-Gawlinski system, which is displayed in nondimensionalized form as follows:
\begin{equation}
\label{eqn:multisystem}
\left\{ \begin{aligned}
\partial_t U &= U (1-U) - d\,U W\\
\partial_t \!V &= r V (1-V) + \alpha \nabla\!\cdot \!\big[ (1-U) \nabla V \big]\\
\partial_t \!W &= c (V-W) + \nabla\!\cdot \!\big( \mathbb{A} \nabla W \big)
\end{aligned} \right.
\end{equation}

for $0\leq t \leq T$ and a suitable domain $\Omega \subset \mathbb{R}^2$, so that the first order differential operator is denoted by $\nabla=(\partial_x,\,\partial_y)$ and $\mathbb{A}$ is an anisotropic, heterogeneous (positive definite) diffusion tensor, namely
\begin{equation}
\label{eqn:anisodiff}
\mathbb{A}(x,y) = \begin{bmatrix}
a(x,y) & b(x,y) \\ b(x,y) & c(x,y)
\end{bmatrix}\,.
\end{equation}


For the sake of convenience, the problem is defined on the square of sides $L$ corresponding to a two-dimensional support typically employed in biology for cell cultures (see Figure~\ref{figure1}(left)), and Neumann boundary conditions are taken into account.

\begin{figure}
\sidecaption
\qquad \includegraphics[scale=.65]{fvca9-simeoni-figure1a}
\qquad \includegraphics[scale=.30]{fvca9-simeoni-figure1b}
\caption{Squared \textit{Petri dish} for in vitro experiments (left) and staggered Cartesian mesh (right)}
\label{figure1}
\end{figure}

\noindent The dynamical variables of system~\eqref{eqn:multisystem} are the scaled unknowns $U(t;x,y)$, $V(t;x,y)$ and $W(t;x,y)$ which represent the normal cells density, the tumor cells density and the extracellular lactic acid concentration in excess, respectively. The cellular populations obey a logistic growth with normalized carrying capacities, and $r$ is the relative growth rate of malignant species. The degradation of healthy tissue by means of the lactic acid is based on the \textit{law of mass action}~\cite{Murr2002} in the form of a second-order kinetic reaction with (adimensionalized) velocity $d$, such a term in the first equation translating the outcome of the pH lowering strategy adopted by cancer cells to enhance their invasive processes.


Another remarkable feature is the density-limited, and possibly degenerate, diffusion of the tumor cells in the second equation, as a consequence of the defense mechanism of confinement by the host tissue experimentally noticed in various contexts~\cite{SwRa2009}, with $\alpha$ being the unconstrained diffusion constant.

Finally, in the third equation, the parameter $c$ encloses both the rate of production for the tumor-associated acidity and a physiological reabsorption rate.


The novelty of system~\eqref{eqn:multisystem} with respect to the original Gatenby-Gawlinski model is that the (excess of) lactic acid diffuses across the microenvironment according to an anisotropic and heterogeneous tensor~\eqref{eqn:anisodiff}, which is suggested by the intrinsic geometrical and/or structural complexity of the extracellular matrix~\cite{Pain2009,DyGr2016}, among other factors. On the other hand, the chemical aggression at the expenses of the local healthy tissue usually includes protease reactions inducing lysis of the extracellular matrix~\cite{ChAn2003}, and therefore the assumption of an isotropic diffusion for cancer cells is also justified.


The environment-driven rate of invasion of malignant species is a widely reported phenomenon~\cite{SwAl2000,ChHi2007,PaHi2013} and we aim at initiating a systematic study of the effect of anisotropic, heterogeneous diffusion operators also within the context of acid-mediated tumor invasion, for which reliable and efficient numerical simulations are necessary because of the unavailability of theoretical results.


\section{Finite volume discretization on staggered grids}
\label{sec:numerics}

The numerical simulation of system~\eqref{eqn:multisystem} relies on finite volume schemes for the spatial operators and the explicit Euler method for time integration, to take advantage from desirable properties like conservativity and high-order accuracy~\cite{LeVe2002}.


In practical applications, since the diffusion tensor~\eqref{eqn:anisodiff} and problem data are derived from the experimental setting, the solution is normally sampled at the nodes of a fixed two-dimensional Cartesian mesh (see Figure~\ref{figure1}(right)). We consider the domain $\Omega=[0,1]\times[0,1]$ which is partitioned into cells of characteristic sizes $\Delta x$ and $\Delta y$, so that the grid points are given by $(x_i,y_j)=(i\,\Delta x,\,j \Delta y)$, $i,j=1,2,...$

\begin{remark}
As already mentioned in~\cite{BeAf1996}, even for computational domains with a complex geometry, using regularly-shaped Cartesian embedded grids generally makes the numerical methods efficient~\cite{BeHe2003,LuFa2016} and indeed easier to implement in parallel compared with boundary fitted structured or unstructured meshes.

\end{remark}

\noindent For the time discretization, starting from $t_0=0$ we proceed by selecting a nonuniform step-size $\Delta t_k$ which is computed after each iteration in order to ensure numerical stability of the finite volume scheme, so that $t_{k+1}=t_k+\Delta t_k$, $k=0,1,...$\\

Finally, we prescribe auxiliary interfacial coordinates, namely $x_{i+\frac12}\!=\!\dfrac{x_i+x_{i+1}}2$ and $y_{j+\frac12}\!=\!\dfrac{y_j+y_{j+1}}2$, to construct dual cells $C_{ij}\!=\!\big(x_{i-\frac12},x_{i+\frac12}\big)\!\times\!\big(y_{j-\frac12},y_{j+\frac12}\big)$ where the numerical unknowns are defined at each time frame as cell-centered approximations
\begin{equation*}

U_{ij}^k \approx \frac1{\big| C_{ij} \big|} \int_{C_{ij}}\hspace{-5pt}U(t_k;x,y)\,dx\,dy\,, \quad \big| C_{ij} \big| = \Delta x\Delta y\,,

\end{equation*}
for instance, and analogously for the other variables of system~\eqref{eqn:multisystem},

together with the piecewise constant projections of the diffusion rates~\eqref{eqn:anisodiff} on the staggered mesh, which are denoted by $a_{ij}$, $b_{ij}$ and $c_{ij}$, for $i,j=1,2,...$


The simplest form of fully discretized problem is presented as follows:
\begin{equation}
\label{eqn:multischeme}
\left\{ \begin{aligned}
U_{ij}^{k+1} &= \,U_{ij}^k + \Delta t_k U_{ij}^k \big(1-U_{ij}^k \big) - d\,\Delta t_k U_{ij}^k W_{ij}^k\\
V_{ij}^{k+1} &= \,V_{ij}^k + r\,\Delta t_k V_{ij}^k \big(1-V_{ij}^k \big) + \alpha\,\Delta t_k \,\mathbb{D}_1(V;U)^k_{ij}\\

W_{ij}^{k+1} &= \,W_{ij}^k + c\,\Delta t_k \big(V_{ij}^k-W_{ij}^k \big) + \Delta t_k \,\mathbb{D}_2(W;\mathbb{A})^k_{ij}

\end{aligned} \right.
\end{equation}
and then this scheme is applied with Neumann-type boundary conditions.

To make the notation more compact, we represent the numerical diffusion operators $\mathbb{D}_1(V;U)^k_{ij}$ and $\mathbb{D}_2(W;\mathbb{A})^k_{ij}$ in~\eqref{eqn:multischeme} by means of the two-dimensional stencils reported in Table~\ref{tab:tumorstencil} and Table~\ref{tab:acidstencil} for the second and third equation, respectively, where the central entries are the coefficients of $V^k_{ij}$ and $W^k_{ij}$, the right-hand side applies to $V^k_{i+1,j}$ and $W^k_{i+1,j}$, the upper side corresponds to $V^k_{i,j+1}$ and $W^k_{i,j+1}$, and similarly for all the others.


\renewcommand{\arraystretch}{2.75}
\begin{table}
\begin{center}
\begin{tabular}{c|c|c}
$0$ & $\dfrac1{\Delta y^2} - \dfrac{U^k_{ij}+U^k_{i,j+1}}{2 \Delta y^2}$ & $0$ \\[3.75pt] \hline 

$\dfrac1{\Delta x^2} - \dfrac{U^k_{i-1,j}+U^k_{ij}}{2 \Delta x^2}\;$ 
& $\; -\dfrac2{\Delta x^2} + \dfrac{U_{i-1,j}+2\,U_{ij}+U_{i+1,j}}{2 \Delta x^2} \qquad$ 
& $\;\dfrac1{\Delta x^2} - \dfrac{U^k_{ij}+U^k_{i+1,j}}{2 \Delta x^2}$ 
\tabularnewline
& $\qquad -\dfrac2{\Delta y^2} + \dfrac{U_{i,j-1}+2\,U_{ij}+U_{i,j+1}}{2 \Delta y^2} \;$ & \\[3.75pt] \hline 


$0$ & $\dfrac1{\Delta y^2} - \dfrac{U^k_{i,j-1}+U^k_{ij}}{2 \Delta y^2}$ & $0$ \\[3.75pt]

\end{tabular}
\end{center}
\caption{Five-point stencil of the standard discretization for the nonlinear tumor cells diffusion.}
\label{tab:tumorstencil}
\end{table}

\renewcommand{\arraystretch}{2.5}
\begin{table}
\begin{center}
\begin{tabular}{c|c|c}
$\; -\dfrac{b_{i-1,j}+b_{i,j+1}}{4 \Delta x \Delta y} \;$ 
& $\dfrac{c_{ij}+c_{i,j+1}}{2 \Delta y^2}$ 
& $\dfrac{b_{i+1,j}+b_{i,j+1}}{4 \Delta x \Delta y}$ \\[5pt] \hline 
$\dfrac{a_{i-1,j}+a_{ij}}{2 \Delta x^2}$ 
& $\; -\dfrac{a_{i-1,j}+2\,a_{ij}+a_{i+1,j}}{2 \Delta x^2}-\dfrac{c_{i,j-1}+2\,c_{ij}+c_{i,j+1}}{2 \Delta y^2} \;$ 
& $\dfrac{a_{ij}+a_{i+1,j}}{2 \Delta x^2}$ \\[5pt] \hline 
$\dfrac{b_{i-1,j}+b_{i,j-1}}{4 \Delta x \Delta y}$ 
& $\dfrac{c_{i,j-1}+c_{ij}}{2 \Delta y^2}$ 
& $\; -\dfrac{b_{i+1,j}+b_{i,j-1}}{4 \Delta x \Delta y} \;$ \\[5pt]
\end{tabular}
\end{center}
\caption{Nine-point stencil of the standard discretization for the (linear) anisotropic, heterogeneous lactic acid dissusion.}
\label{tab:acidstencil}
\end{table}

\noindent The consistency and stability of finite volume schemes on staggered Cartesian grids for anisotropic and heterogeneous diffusion tensors is extensively treated in~\cite{DaSi2014}.

As a matter of fact, these questions are closely related to the discrete maximum principle for the \textit{Beltrami color flow}, which is typically used in image processing~\cite{DaDi2007}.\\


Although the major drawback of explicit algorithms for reaction-diffusion systems like~\eqref{eqn:multischeme} resides precisely in the stringent limitations on the time-step enforced by both the mesh geometry and diffusion rates, through the \textit{CFL-condition}, the parallel implementation using GPUs rehabilitates such methods for the benefits obtained in terms of drastic reduction of the computational time, in contrast to the expensive simulations produced with sequential codes (refer to Section~\ref{sec:performance}).


Moreover, the accuracy in designing spatially compact numerical operators is motivated by the requirements of an optimal parallel implementation~\cite{Fels2008,KlWa2009}, since the nearest-neighbor communication standard is extremely fast with the need for small amounts of local storage in the sub-processors, because only few values are involved to update the numerical solution at each grid node (as shown in Table~\ref{tab:tumorstencil} and Table~\ref{tab:acidstencil}).


\begin{remark}
We point out that the extra-diagonal entries in Table~\ref{tab:acidstencil}, which arise from the anisotropic structure of diffusion tensors~\eqref{eqn:anisodiff}, are indeed responsible for an increased computational complexity, since the matrix-valued version of algorithm~\eqref{eqn:multischeme} does not benefit from sparsity-inducing properties. On the other hand, a different approach to efficient computational arrangements is required for the five-point stencil in Table~\ref{tab:tumorstencil}, because the corresponding matrix-valued scheme is actually nonlinear and need to be reconstructed at each time iteration. Nevertheless, various modern C++ libraries for numerically solving ordinary differential equations are presently available, which are compatible with running on CUDA GPUs through the Thrust interface -- \url{http://thrust.github.io}

\end{remark}

A trademark of the Gatenby-Gawlinski model consists in the capability of reproducing the biologically relevant phenomenon of propagating fronts and, in particular, to predict the development of a \textit{hypocellular interstitial gap}, that is a local region practically depleted of cells where both the healthy and cancerous tissue densities are negligible, which has been experimentally observed~\cite{GaGa1996}.


In the context of acid-mediated tumor invasion, the dynamical configurations of system~\eqref{eqn:multisystem} with respect to the reaction parameter $d$ essentially exhibit two types of qualitative behavior:
\begin{enumerate}

\item{the \textit{heterogeneous invasion}, occurring for moderate chemical aggression ($d<1$) and characterized by the coexistence of malignant and normal tissues behind the propagating tumor front;}


\item{the more aggressive \textit{homogeneous invasion}, which refers to the complete degradation ($d>1$) of local healthy tissue by means of the lactic acid ahead of the propagating tumor front (see Figure~\ref{figure2}(left)); in that case, the emergence of an interstitial gap and its geometrical structure are submitted also to the magnitude of the diffusion rates~\eqref{eqn:anisodiff}, and larger values of $d$ are needed to compensate for lower acid diffusivity in some directions (see Figure~\ref{figure2}(right)).}

\end{enumerate}


\begin{figure}
\sidecaption
\includegraphics[scale=.165]{fvca9-simeoni-figure2a}
\qquad \includegraphics[scale=.165]{fvca9-simeoni-figure2b}
\caption{Interstitial gap between cancerous and healthy cells for lower (left) and higher (right) values of the aggressiveness parameter $d$ (in the case of lower acid diffusivity along the vertical direction)}
\label{figure2}
\end{figure}


These numerical results suggest that anisotropic and heterogenous acid diffusion is determinant in shaping the propagating tumor fronts, and further investigations are in order for recovering experimental observations about the acid-mediated invasive processes.


In this respect, finite volume schemes on staggered grids are particularly well-suited for their ability in accurately reproducing multi-dimensional interfacial problems, especially together with GPU-based computing for biological applications~\cite{MaMi2013,PeMa2019}, and therefore comparable to finite element simulations~\cite{MoSi2019}.


\section{GPU programming and performance evaluation}
\label{sec:performance}

Beside possibly anomalous regularizing effects, a crucial feature for the mathematical modeling of invasion processes with anisotropic and heterogeneous diffusion tensors~\eqref{eqn:anisodiff} is the emergence of strikingly nontrivial patterns~\cite{HiPa2013}.


Consequently, the numerical simulations often requires high spatial resolution, and then long-lasting computing execution, for the large amount of data to be traded in order to accurately capture the details of biological phenomena. Moreover, especially for clinical practitioners and applied scientists involved in setting up realistic experiments, the possibility of running fast comparative simulations using simple algorithms implemented into affordable processors is of primary interest.


Parallel computing based on modern Graphics Processing Units (GPUs) has the advantage of high performance at a relatively low energetically and monetary costs.\\


The codes used in this article are programmed using the NVIDIA$^{\textcopyright}$ Compute Unified Device Architecture (CUDA) platform, which is designed to support GPUs execution for data parallelization~\cite{CUDA2019}. Through the CUDA implementation, graphics cards are programmed with a medium-level language, which is recognized as an extension to C/C++, without requiring sophisticated hardware expertise~\cite{KiHw2016,SaKa2010}.\\

All simulations\footnote{The code for reproducing the numerical tests is available upon request to the authors.} in this article are performed using NVIDIA$^{\textcopyright}$ graphics cards GTX 1080 with 2560 CUDA cores and 8 Gb RAM, installed on a processor HP DL585G7 4 AMD Opteron 6128 with 8 cores, clock frequency 2.0 GHz, 64 Gb RAM, operating system Linux centOS 6.5 amd64, compiler GNU gcc 4.4.7 and NVIDIA$^{\tiny{\textcopyright}}$ CUDA 9.1 Linux 64 bit toolbox. For a comparative study of computing performance, the same algorithms are implemented serially on a single CPU processor HP DL585G7 4 AMD Opteron 6128 CPU with 8 cores.


The standard logical steps for implementing the parallel codes are as follows:
\begin{enumerate}

\item{the initial data and numerical parameters are loaded into the CPU (host) memory, and then transferred to the GPU (device) storage, namely the \textit{global memory}};
\item{at each time iteration, the GPU provides massive computing activities, by executing in parallel across the sub-processors the algorithm for the spatial discretization, since the current unknowns only depend on those already computed at the previous iteration};
\item{the GPU uploads the local simulation results to the CPU memory};
\item{the CPU updates the value of the time-step, according to the constraints imposed to guarantee the numerical stability, and then it restarts/stops the parallel computing process.}

\end{enumerate}

\noindent From the point of view of programming, this approach leads to designate an external \textit{loop} for the time integration, whereas the spatial approximation at each iteration is carried out in parallel by the GPU within a \textit{master-slave} model, which is sketched in Figure~\ref{figure3}(left). Finally, the performance comparison between serial and parallel execution for the above simulations is shown in Figure~\ref{figure3}(right), where the \textit{speed-up} is illustrated as a function of the mesh size, for instance.


\begin{figure}
\sidecaption
\includegraphics[scale=.30]{fvca9-simeoni-figure3a} 
\quad \includegraphics[scale=.45]{fvca9-simeoni-figure3b}
\caption{General steps of GPU-based parallel codes (left) and \textit{speed-up} evaluation (right)}
\label{figure3}
\end{figure}


The CUDA implementation of the numerical codes is currently undertaking further optimizations, for example by employing other GPU \textit{local memories} for achieving \textit{dynamic parallelism}, and these issues are oriented to the three-dimensional modeling of experimental data.


\begin{acknowledgement}

This work is supported by the French government, managed by the ANR under the UCA JEDI Investments for the Future project, reference n. ANR-15-IDEX-01. All numerical simulations are performed on the Linux HPC parallel cluster Caliban -- \url{https://caliband.disim.univaq.it} -- at the University of L'Aquila (Italy). This article is greatly motivated by fruitful discussions with the Systems Biology Group -- \url{https://sbglab.org} -- at the Department of Experimental Medicine, Sapienza University of Rome (Italy).


\end{acknowledgement}

\input{fvca9-simeoni-references}

\end{document}
{"config": "arxiv", "file": "2006.01748/fvca9-simeoni.tex"}