A hybrid discrete-continuum modelling approach for the interactions of the immune system with oncolytic viral infections

TL;DR

The results highlight that a too rapid immune response, before the infection is well-established, appears to decrease the efficacy of the therapy and thus some care is needed when oncolytic virotherapy is combined with immunotherapy, suggesting the importance of clinically improving the modulation of the immune response according to the tumour's characteristics and to the immune capabilities of the patients.

Abstract

Oncolytic virotherapy, utilizing genetically modified viruses to combat cancer and trigger anti-cancer immune responses, has garnered significant attention in recent years. In our previous work arXiv:2305.12386, we developed a stochastic agent-based model elucidating the spatial dynamics of infected and uninfected cells within solid tumours. Building upon this foundation, we present a novel stochastic agent-based model to describe the intricate interplay between the virus and the immune system; the agents' dynamics are coupled with a balance equation for the concentration of the chemoattractant that guides the movement of immune cells. We formally derive the continuum limit of the model and carry out a systematic quantitative comparison between this system of PDEs and the individual-based model in two spatial dimensions. Furthermore, we describe the traveling waves of the three populations, with the uninfected proliferative cells trying to escape from the infected cells while immune cells infiltrate the tumour. Simulations show a good agreement between agent-based approaches and numerical results for the continuum model. Some parameter ranges give rise to oscillations of cell number in both models, in line with the behaviour of the corresponding nonspatial model, which presents Hopf bifurcations. Nevertheless, in some situations the behaviours of the two models may differ significantly, suggesting that stochasticity plays a key role in the dynamics. Our results highlight that a too rapid immune response, before the infection is well-established, appears to decrease the efficacy of the therapy and thus some care is needed when oncolytic virotherapy is combined with immunotherapy. This further suggests the importance of clinically improving the modulation of the immune response according to the tumour's characteristics and to the immune capabilities of the patients.

Figures (12)

• Figure 1: Schematic representation of the rules governing cell dynamics in the stochastic models. Uninfected cells are represented in blue, infected cells in red and immune cells in green. Uninfected cells may proliferate or die according to the total density, move, become infected upon contact with infected cells and die upon contact with immune cells. Infected cells may move, die with constant probability and die upon contact with immune cells. Immune cells may enter the domain, move with the probabilities given in Eq. \ref{['eq:F']} and die with constant probability. The model also considers the dynamics of the chemoattractant, which are not included in the figure due to the different modelling approach adopted (i.e., density-based and deterministic instead of individual-based and stochastic).
• Figure 2: One parameter bifurcations in $\alpha$, $\zeta$ and $\beta$ of Eq. \ref{['eq:ode']}, with other parameters as in Table \ref{['tab:parameters']}. The immune killing rate $\zeta$ has been set to the base value $0.50\;$h$^{-1}$. In order to facilitate comparison with the forthcoming two-dimensional simulations, we set $\alpha=\pi r^2 \alpha_z$ with $r=5\;$mm (corresponding to a late stage of tumour growth). The green dots show the maximum and minimum values of $u$ during the oscillations of the stable limit cycle. The solid lines show the value of the equilibrium of $u$; the line is red if the equilibrium is stable and black if it is unstable. Hopf bifurcations are denoted by HB. Observe that for low values of $\beta$ the infection-free equilibrium close to carrying capacity is stable.
• Figure 3: Numerical simulation of Eq. \ref{['eq:ode']} with the parameters as in Table \ref{['tab:parameters']} and different values of the immune killing rate $\zeta$. As in Fig. \ref{['fig:bif_onepar']}, we set $\alpha=\pi r^2 \alpha_z$ with $r=5\;$mm. Uninfected tumour cells are plotted in blue, infected tumour cells in red and immune cells in green. The oscillations become wider as $\zeta$ increases, in accordance with the bifurcation diagram of Fig. \ref{['fig:bif_onepar']}b.
• Figure 4: Comparisons of the scenarios with weak immune response and no infection (solid lines, $\zeta=0.50\,$h$^{-1}$ and $R_i=0$), strong immune response and no infection (dashed lines, $\zeta=5.00\,$h$^{-1}$ and $R_i=0$), weak immune response and central infection (dash-dotted lines, $\zeta=0.50\,$h$^{-1}$). All the other parameters take the values given in Table \ref{['tab:parameters']}. The results of the agent-based model are averaged over ten simulations. Panels (a-b) represent the horizontal section of the domain $[-L,L]\times\{0\}$: uninfected cells are in blue, infected cells in red (only shown with dash-dotted lines, as in the other cases there is no infection) and chemokines in yellow. The vertical blue dashed lines represent the expected positions of the uninfected invading front, travelling at speed $2\sqrt{D_u p}$. Panel (c) shows the sum of tumour cells over time. The results of the agent-based model in the three scenarios are represented with solid, dashed and dash-dotted lines (as explained in the legend). The results of the continuum models are represented with dotted lines in all the three cases, as there is an excellent agreement with the stochastic counterparts.
• Figure 5: A single numerical simulation of the agent-based model with the parameters given in Table \ref{['tab:parameters']}, $\zeta=0.50\,$h$^{-1}$ and wide oncolytic viral infection (i.e., $R_i=R_u$). The dashed cyan circles in panels (a) and (c) represent the expected positions of the tumour invading front, travelling at speed $2\sqrt{D_u p}$. The dotted green circles in panels (a), (b) and (c) represent the internal minimum of the numerical solution of Eq. \ref{['eq:pderad']}. The dashed red circle in panel (b) represents the front given by the numerical solution of Eq. \ref{['eq:pderad']}. In panel (d), solid lines refer to the agent-based model (uninfected, infected and immune cells are represented respectively in blue, red and green) and dotted lines refer to the continuum model. In all the cases the maximum of the axes and the colorbars correspond to the maximum over time of the quantity plotted.