<figure><img src="image_1.png" /></figure>

HAL
open science

# A continuous approach of modeling tumorigenesis and axons regulation for the pancreatic cancer

Marie-Jose Chaaya, Sophie Chauvet, Florence Hubert, Fanny Mann, Mathieu Mezache, Pierre Pudlo

<figure><img src="image_2.png" /></figure>

## ► To cite this version:

Marie-Jose Chaaya, Sophie Chauvet, Florence Hubert, Fanny Mann, Mathieu Mezache, et al. A continuous approach of modeling tumorigenesis and axons regulation for the pancreatic cancer. 2024. [hal-04527648](https://hal.inrae.fr/hal-04527648)

HAL Id: hal-04527648

https://hal.inrae.fr/hal-04527648

Preprint submitted on 2 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# A continuous approach of modeling tumorigenesis and axons regulation for the pancreatic cancer

Marie-Jose Chaaya, Sophie Chauvet, Florence Hubert, Fanny Mann, Mathieu Mezache, Pierre Pudlo

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Marie-Jose Chaaya, Sophie Chauvet, Florence Hubert, Fanny Mann, Mathieu Mezache, et al. A continuous approach of modeling tumorigenesis and axons regulation for the pancreatic cancer. Journal of Theoretical Biology, 2024, 595, pp.111967. 10.1016/j.jtbi.2024.111967. hal-04527648v2

HAL Id: hal-04527648

https://hal.inrae.fr/hal-04527648v2

Submitted on 5 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A continuous approach of modeling tumorigenesis and axons regulation for the pancreatic cancer.

Marie-Jose Chaaya<fnref n="2" />, Sophie Chauvet<fnref n="1" />, Florence Hubert<fnref n="2" />, Fanny Mann<fnref n="1" />, Mathieu Mezache<fnref n="2" />,<fnref n="3" />, and Pierre Pudlo<fnref n="2" />

**Keywords:** Dynamical system, Partial differential equations, Cancer, Parameter calibration, In silico denervation

## Abstract

The pancreatic innervation undergoes dynamic remodeling during the development of pancreatic ductal adenocarcinoma (PDAC). Denervation experiments have shown that different types of axons can exert either pro- or anti-tumor effects, but conflicting results exist in the literature, leaving the overall influence of the nervous system on PDAC incompletely understood. To address this gap, we propose a continuous mathematical model of nerve-tumor interactions that allows in silico simulation of denervation at different phases of tumor development. This model takes into account the pro- or anti-tumor properties of different types of axons (sympathetic or sensory) and their distinct remodeling dynamics during PDAC development. We observe a "shift effect" where an initial pro-tumor effect of sympathetic axon denervation is later outweighed by the anti-tumor effect of sensory axon denervation, leading to a transition from an overall protective to a deleterious role of the nervous system on PDAC tumorigenesis. Our model also highlights the importance of the impact of sympathetic axon remodeling dynamics on tumor progression. These findings may guide strategies targeting the nervous system to improve PDAC treatment.

## 1 Introduction

The nervous system plays an important role in regulating various bodily functions and disease processes, including cancer development and progression [Winkler et al., 2023]. Denervation studies have shown that the peripheral nervous system (PNS) can either promote or inhibit cancer growth, depending on the specific types of nerves and cancers involved. For example, in mouse models of pancreatic ductal adenocarcinoma (PDAC), selective ablation of pancreatic sympathetic innervation has been associated with accelerated tumor growth, increased metastasis and decreased survival [Guillot et al., 2022]. Conversely, removal or silencing of sensory neurons has been shown to slow tumor growth and improve survival ( [Saloman et al., 2016]; [Sinha et al., 2017]). These findings, together with others, suggest a broad model in which the autonomic nervous system (including both its sympathetic and parasympathetic branches) has an anti-tumor property in PDAC, whereas the sensory nervous system has an inverse pro-tumor activity (Figure 1).

According to the data presented, the combined effect of the PNS on PDAC tumour progression results from a mixture of pro- and anti-tumour activities. This combined effect can be assessed in experimental models by surgical denervation of the pancreas, which disrupts the mixed nerve of sympathetic and sensory axons supplying the pancreas. Such an intervention mimics procedures performed in patients (eg. celiac neurolysis and splanchnicectomy) to manage abdominal pain in unresectable PDAC. However, surgical denervation studies in mice have led to divergent results. A pro-tumor effect was observed when denervation was performed before the onset of pathology [Guillot et al., 2022], whereas an anti-tumor effect was observed when denervation was performed after tumor establishment [Renz et al., 2018a]. The reasons behind these conflicting results are not yet understood. One hypothesis suggests a switch in sympathetic function over the course of pathology. Initially, sympathetic axons may have an anti-tumor effect in the early (pre-)cancer stages, before possibly switching to a promoting role on tumor growth at a later stage. Another interpretation is that the relative abundance of pro- and anti-tumoral axons may vary at different denervation times. Indeed, studies have shown different remodeling patterns for sympathetic and sensory innervation. Higher levels of sympathetic innervation are found in pre-cancerous lesions compared to both healthy and cancerous tissues, while sensory axons are markedly increased within the cancerous

<footnote n="1">¹Aix Marseille Univ, CNRS, IBDM (UMR 7288), Turing Centre for Living systems, Marseille, France</footnote>

<footnote n="2">²Aix Marseille Univ, CNRS, I2M (UMR 7373), Turing Centre for Living systems, Marseille, France</footnote>

<footnote n="3">³Université Paris-Saclay, INRAE, MAIAGE (UR 1404), 78350 Jouy-en-Josas, France</footnote>
<figure><img src="image_3.png" /></figure>

# A continuous approach of modeling tumorigenesis and axons regulation for the pancreatic cancer

Marie-Jose Chaaya<sup>b</sup>, Sophie Chauvet<sup>a</sup>, Florence Hubert<sup>b</sup>, Fanny Mann<sup>a</sup>, Mathieu Mezache<sup>b,c,*</sup>, Pierre Pudlo<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Dynamical system

Partial differential equations

Cancer

Parameter calibration

In silico denervation

## ABSTRACT

The pancreatic innervation undergoes dynamic remodeling during the development of pancreatic ductal adenocarcinoma (PDAC). Denervation experiments have shown that different types of axons can exert either pro- or anti-tumor effects, but conflicting results exist in the literature, leaving the overall influence of the nervous system on PDAC incompletely understood. To address this gap, we propose a continuous mathematical model of nerve-tumor interactions that allows in silico simulation of denervation at different phases of tumor development. This model takes into account the pro- or anti-tumor properties of different types of axons (sympathetic or sensory) and their distinct remodeling dynamics during PDAC development. We observe a "shift effect" where an initial pro-tumor effect of sympathetic axon denervation is later outweighed by the anti-tumor effect of sensory axon denervation, leading to a transition from an overall protective to a deleterious role of the nervous system on PDAC tumorigenesis. Our model also highlights the importance of the impact of sympathetic axon remodeling dynamics on tumor progression. These findings may guide strategies targeting the nervous system to improve PDAC treatment.

## 1. Introduction

The nervous system plays an important role in regulating various bodily functions and disease processes, including cancer development and progression (Winkler et al., 2023). Denervation studies have shown that the peripheral nervous system (PNS) can either promote or inhibit cancer growth, depending on the specific types of nerves and cancers involved. For example, in mouse models of pancreatic ductal adenocarcinoma (PDAC), selective ablation of pancreatic sympathetic innervation has been associated with accelerated tumor growth, increased metastasis and decreased survival (Guillot et al., 2022). Conversely, removal or silencing of sensory neurons has been shown to slow tumor growth and improve survival (Saloman et al., 2016; Sinha et al., 2017; Liddle, 2007). These findings, together with others, suggest a broad model in which the autonomic nervous system (including both its sympathetic and parasympathetic branches) has an anti-tumor property in PDAC, whereas the sensory nervous system has an inverse pro-tumor activity (Fig. 1).

According to the data presented, the combined effect of the PNS on PDAC tumor progression results from a mixture of pro- and anti-tumor activities. This combined effect can be assessed in experimental

models by surgical denervation of the pancreas, which disrupts the mixed nerve of sympathetic and sensory axons supplying the pancreas. Such an intervention mimics procedures performed in patients (eg. celiac neurolysis and splanchnicectomy) to manage abdominal pain in unresectable PDAC. However, surgical denervation studies in mice have led to divergent results. A pro-tumor effect was observed when denervation was performed before the onset of pathology (Guillot et al., 2022), whereas an anti-tumor effect was observed when denervation was performed after tumor establishment (Renz et al., 2018a). The reasons behind these conflicting results are not yet understood. One hypothesis suggests a switch in sympathetic function over the course of pathology. Initially, sympathetic axons may have an anti-tumor effect in the early (pre-)cancer stages, before possibly switching to a promoting role on tumor growth at a later stage. Another interpretation is that the relative abundance of pro- and anti-tumoral axons may vary at different denervation times. Indeed, studies have shown different remodeling patterns for sympathetic and sensory innervation. Higher levels of sympathetic innervation are found in pre-cancerous lesions compared to both healthy and cancerous tissues, while sensory axons

<footnote><sup>a</sup> Aix Marseille Univ, CNRS, IBDM (UMR 7288), Turing Centre for Living systems, Marseille, France</footnote>

<footnote><sup>b</sup> Aix Marseille Univ, CNRS, I2M (UMR 7373), Turing Centre for Living systems, Marseille, France</footnote>

<footnote><sup>c</sup> Université Paris-Saclay, INRAE, MaIAGE (UR 1404), 78350 Jouy-en-Josas, France</footnote>

<footnote>* Corresponding author at: Université Paris-Saclay, INRAE, MaIAGE (UR 1404), 78350 Jouy-en-Josas, France.
E-mail address: mathieu.mezache@inrae.fr (M. Mezache).</footnote>
lesions ([Guillot et al., 2022]; [Demir et al., 2015]).

These findings highlight the importance of the PNS as a potential therapeutic target for the modulation of PDAC. However, they also highlight the need for a deeper understanding of the individual and combined effects of the different axon types, taking into account their remodeling patterns and influence on the tumor, in order to develop strategies to deplete or inhibit the PNS in PDAC. To address this need, we develop a mathematical model that allow us to study the consequences of denervating sensory and/or sympathetic axons at different times during tumor development and progression.

So far, two mathematical models have been developed to better understand the influence of the PNS on cancer. The first model investigated the pro-tumoral effect of the autonomic nervous system in prostate cancer [Lolas et al., 2016]. However, due to the opposite function of the sympathetic innervation in pancreatic cancer, another mathematical model was developed to specifically study PDAC [Chauvet et al., 2023]. This model formalized the interactions between cancer progression and axons using a compartmentalized differential equations model, where each compartment corresponds to a stage of cancer progression (healthy, pre-cancerous, cancerous). The asymptotic behavior of the system shows that the pathological state, where only cancer cells persist, is globally asymptotically stable. The impact of denervation was simulated in silico and recapitulated the biological data of denervation performed at early stage, before the onset of pathology. However, other times of denervation have not been investigated.

The aim of this paper is to introduce a new continuous mathematical model describing the relationship between axons and cancer, and to simulate in silico denervation of sympathetic and sensory neurons at different times. The continuous model considers the cell phenotype (from healthy to cancerous) as a continuous variable. Consequently, the description of cancer progression and axon remodeling occurring during this process becomes more precise.

We organized the paper as follows. In Section 2.1, we set up the model illustrating the effect that PNS axons have on cancer development and progression. A theoretical study on the well-posedness of the model is given in Section 2.2. This is followed by a qualitative study in Section 2.3 where we show that the final state of the disease can belong to one out of three cases: either no cancer cells exist, either cancer and healthy cells co-exist, or only cancer cells exist. Also, explicit bounds on the time of appearance of the first cancer cells are given under some hypotheses on the parameters. Next, in Section 3, we give precision on the parameters of the model and explain how to denervate in silico. We then calibrate the model in Section 4 by minimizing a 2-dimensional criterion. We obtain several sets of parameters that fit the data. Using some of the sets of parameters obtained, we apply numerical simulations to study in silico denervation using an indicator of invasive potential detailed in Section 5.1. In Section 5.2, we perform in silico denervation separately for each axon type (sympathetic or sensory) or simultaneously at defined time points. The dynamical system and its in silico denervation are implemented using a finite volume approach and the algorithm can be found at (https://github.com/MarieJosec/PDE_Axons_Innerv). In Section 5.3, denervation is performed for each axon type or for both at different time points. The results recapitulate the different outcomes of surgical denervation observed in [Guillot et al., 2022] and [Renz et al., 2018a] [Renz et al., 2018a] and support a 'shift effect' of PNS function in PDAC from anti-tumor to pro-tumor. The model predicts that this transition does not occur through a change in sympathetic function, but rather depends on the remodeling dynamics of the sympathetic and sensory axons and on the strength of the sympathetic inhibition on tumor growth.

# 2 Mathematical model

In this section, we establish and explore an original continuous model for studying cancer progression and its regulation by axons. Initially, a compartmental model was introduced in [Chauvet et al., 2023] to investigate the interactions. However, this model, based on ordinary differential equations (ODEs), focused on discrete stages of cancer progression from a healthy state to a tumoral stage. Each stage was treated as a variable, compartmentalizing the different steps, which failed to capture the continuous nature of the biological process. To address this limitation, we introduce another mathematical model that adopts a continuous framework for cell phenotypes, spanning from a healthy state to a tumoral stage, thereby bringing it closer to the reality of the biological process. This model also considers the heterogeneity of the cell phenotype. In alignment with recent literature (cf. [Eftimie and Gibelli, 2020], which introduces a phenotype-structured model to describe the heterogeneous population of macrophages), we propose employing structured population dynamics in the model.

## 2.1 Model presentation

We focus on the pancreas as the domain of our model, by taking into account the pancreatic cells and the nerve axons. In this model, we denote the unknowns $Q, A_1$ and $A_2$ where:
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Graphical description of the state-of-the-art concerning the interactions between sensory axons, sympathetic axons and tumor cells (PDAC). The arrow with the symbol → denotes a promoting effect and the arrow with the symbol ← denotes an inhibiting effect.</figcaption></figure>

are markedly increased within the cancerous lesions (Guillot et al., 2022; Demir et al., 2015).

These findings highlight the importance of the PNS as a potential therapeutic target for the modulation of PDAC. However, they also highlight the need for a deeper understanding of the individual and combined effects of the different axon types, taking into account their remodeling patterns and influence on the tumor, in order to develop strategies to deplete or inhibit the PNS in PDAC. To address this need, we develop a mathematical model that allow us to study the consequences of denervating sensory and/or sympathetic axons at different times during tumor development and progression.

So far, two mathematical models have been developed to better understand the influence of the PNS on cancer. The first model investigated the pro-tumoral effect of the autonomic nervous system in prostate cancer (Lolas et al., 2016). However, due to the opposite function of the sympathetic innervation in pancreatic cancer, another mathematical model was developed to specifically study PDAC (Chauvet et al., 2023). This model formalized the interactions between cancer progression and axons using a compartmentalized differential equations model, where each compartment corresponds to a stage of cancer progression (healthy, pre-cancerous, cancerous). The asymptotic behavior of the system shows that the pathological state, where only cancer cells persist, is globally asymptotically stable. The impact of denervation was simulated in silico and recapitulated the biological data of denervation performed at early stage, before the onset of pathology. However, other times of denervation have not been investigated.

The aim of this paper is to introduce a new continuous mathematical model describing the relationship between axons and cancer, and to simulate in silico denervation of sympathetic and sensory neurons at different times. The continuous model considers the cell phenotype (from healthy to cancerous) as a continuous variable. Consequently, the description of cancer progression and axon remodeling occurring during this process becomes more precise.

We organized the paper as follows. In Section 2.1, we set up the model illustrating the effect that PNS axons have on cancer development and progression. A theoretical study on the well-posedness of the model is given in Section 2.2. This is followed by a qualitative study in Section 2.3 where we show that the final state of the disease can belong

to one out of three cases: either no cancer cells exist, either cancer and healthy cells co-exist, or only cancer cells exist. Also, explicit bounds on the time of appearance of the first cancer cells are given under some hypotheses on the parameters. Next, in Section 3, we give precision on the parameters of the model and explain how to denervate in silico. We then calibrate the model in Section 4 by minimizing a 2-dimensional criterion. We obtain several sets of parameters that fit the data. Using some of the sets of parameters obtained, we apply numerical simulations to study in silico denervation using an indicator of invasive potential detailed in Section 5.1. In Section 5.2, we perform in silico denervation separately for each axon type (sympathetic or sensory) or simultaneously at defined time points. The dynamical system and its in silico denervation are implemented using a finite volume approach and the algorithm can be found at (https://github.com/MarieJosec/PDE_Axons_Innerv). In Section 5.3, denervation is performed for each axon type or for both at different time points. The results recapitulate the different outcomes of surgical denervation observed in Guillot et al. (2022) and Renz et al. (2018a) [Renz et al. 2018a] and support a 'shift effect' of PNS function in PDAC from anti-tumor to pro-tumor. The model predicts that this transition does not occur through a change in sympathetic function, but rather depends on the remodeling dynamics of the sympathetic and sensory axons and on the strength of the sympathetic inhibition on tumor growth.

## 2. Mathematical model

In this section, we establish and explore an original continuous model for studying cancer progression and its regulation by axons. Initially, a compartmental model was introduced in Chauvet et al. (2023) to investigate the interactions. However, this model, based on ordinary differential equations (ODEs), focused on discrete stages of cancer progression from a healthy state to a tumoral stage. Each stage was treated as a variable, compartmentalizing the different steps, which failed to capture the continuous nature of the biological process. To address this limitation, we introduce another mathematical model that adopts a continuous framework for cell phenotypes, spanning from a healthy state to a tumoral stage, thereby bringing it closer to the reality of the biological process. This model also considers the heterogeneity
<figure><img src="image_3.png" /><figcaption>Figure 1: Graphical description of the state-of-the-art concerning the interactions between sensory axons, parasympathetic neurons, sympathetic axons and tumor cells (PDAC). The arrow with the symbol → denotes a promoting effect and the arrow with the symbol ← denotes an inhibiting effect.</figcaption></figure>

* $Q(t, x)$ is the density of cells structured by the progression state of the disease. The variable $t$ corresponds to time with $t \in [0, +\infty)$ and the variable $x$ corresponds to the phenotype of the cell where $x \in \Omega$ with $\Omega := \mathbb{R}$ or $\Omega := (-L, L)$ with the constant $L > 0$ (finite truncated domain). Cells with a positive phenotype correspond to cancerous cells. Whereas for $x \in (-L, 0)$, we start from healthy cells undergoing precancerous phenotypes before reaching the cancerous state. The lower the phenotype variable $x$ is, the healthier the cells $Q(t, x)$ are.

* The variable $A_1(t)$ is the normalized density of the sympathetic axons with respect to time $t \in [0, +\infty)$.

* The variable $A_2(t)$ is the normalized density of sensory axons with respect to time $t \in [0, +\infty)$.

Hence, the variables $A_i$ for $i = 1, 2$ are unit-less and are bounded, i.e. $A_i(t) \in (0, 1)$ for $t \ge 0$.

The total amount of pancreatic cells $N(t)$ and the total amount of cancer cells $N_c(t)$ can be obtained from the density $Q$ by

$$N(t) = \int_{\Omega} Q(t, x)dx, \quad N_c(t) = \int_{\Omega} \psi(x)Q(t, x)dx,$$

where $\psi$ is an indicator function whose support is a subset of $\mathbb{R}_+$. These two macroscopic quantities play significant roles in the dynamics of the model. They allow to compute crucial indicators such as the proportion of cancer cells in the system or the growth of the size of the pancreas induced by cancer cells. For instance, the ratio $N_c(t)/N(t)$ corresponds to the proportion of cancer cells in the system. If the ratio is equal to 1 then all the cells are cancerous ones and if the ratio is equal to 0, then no cancer cells are present in the system. Also the ratio $N(t)/N(0)$ gives information on the growth of the size of the pancreas.

Moreover, we introduce the notation $\mathcal{X} := (Q, A_1, A_2)$ to group the unknowns into a tuple. Hence, the dynamic of $Q(t, x)$ is given by the following transport-growth equation:

$$\partial_t Q(t, x) + \partial_x [f(t, x; \mathcal{X})Q(t, x)] = g(t, x; \mathcal{X}).$$

**Progression of the disease.** We model the evolution of the disease as a transport term on the phenotype axis for the cell densities $Q$ in the partial differential equation governing the dynamics of cells. The speed of progression of the disease (i.e. the transport speed) denoted $f(t, x; \mathcal{X})$ is regulated by the presence of the axons and the cancer cells. It takes the following form

$$f(t, x; \mathcal{X}) := \pi(x) [1 - \beta(x)\rho(A_1(t)) + \delta(x)A_2(t)] + \eta \left(x, \frac{N_c(t)}{N(t)}\right) \quad (1)$$

where
of the cell phenotype. In alignment with recent literature (cf. Eftimie and Gibelli, 2020, which introduces a phenotype-structured model to describe the heterogeneous population of macrophages), we propose employing structured population dynamics in the model.

## 2.1. Model presentation

We focus on the pancreas as the domain of our model, by taking into account the pancreatic cells and the nerve axons. Despite the fact that parasympathetic nerves have a proven effect on pancreatic cancer, cf. Renz et al. (2018b), to our knowledge there are no quantitative data on parasympathetic remodeling in PDAC tumorigenesis. Hence, we choose to focus on modeling sensory/sympathetic versus PDAC interactions (cf. Fig. 1). In this model, we denote the unknowns $Q, A_1$ and $A_2$ where:

* $Q(t, x)$ is the density of cells structured by the progression state of the disease. The variable $t$ corresponds to time with $t \in [0, +\infty)$ and the variable $x$ corresponds to the phenotype of the cell where $x \in \Omega$ with $\Omega := \mathbb{R}$ or $\Omega := (-L, L)$ with the constant $L > 0$ (finite truncated domain defined for the purposes of numerical simulations). Cells with a positive phenotype correspond to cancerous cells. Whereas for $x \le 0$, we start from healthy cells undergoing precancerous phenotypes before reaching the cancerous state. The lower the phenotype variable $x$ is, the healthier the cells $Q(t, x)$ are.

* The variable $A_1(t)$ is the normalized density of the sympathetic axons with respect to time $t \in [0, +\infty)$. We denote by “normalized density” an original density value which has been scaled by dividing it by the maximum density in order to obtain a value on standardized scale.

* The variable $A_2(t)$ is the normalized density of sensory axons with respect to time $t \in [0, +\infty)$.

Hence, the variables $A_i$ for $i = 1, 2$ are unit-less and are bounded, i.e. $A_i(t) \in (0, 1)$ for $t \ge 0$.

The total amount of pancreatic cells $N(t)$ and the total amount of cancer cells $N_c(t)$ can be obtained from the density $Q$ by

$$N(t) = \int_{\Omega} Q(t, x)dx, \quad N_c(t) = \int_{\Omega} \psi(x)Q(t, x)dx,$$

where $\psi$ is an indicator function whose support is a subset of $\mathbb{R}_+$. These two macroscopic quantities play significant roles in the dynamics of the model. They allow to compute crucial indicators such as the proportion of cancer cells or the growth of the size of the pancreas induced by cancer cells. For instance, the ratio $N_c(t)/N(t)$ corresponds to the total amount of cancer cells at time $t$ over the total amount of cells at time $t$. If the ratio is equal to 1 then all the cells are cancerous ones and if the ratio is equal to 0, then no cancer cells are present at time $t$. Also the ratio $N(t)/N(0)$ gives information on the growth of the size of the pancreas.

Moreover, we introduce the notation $\mathcal{X} := (Q, A_1, A_2)$ to group the unknowns into a tuple. Hence, the dynamic of $Q(t, x)$ is given by the following transport-growth equation:

$$\partial_t Q(t, x) + \underbrace{\partial_x [f(t, x; \mathcal{X})Q(t, x)]}_{\text{Progression of the disease}} = \underbrace{g(t, x; \mathcal{X})}_{\text{Growth term}}.$$

**Progression of the disease.** We model the evolution of the disease as a transport term on the phenotype axis for the cell densities $Q$ in the partial differential equation governing the dynamics of cells. The speed of progression of the disease (i.e. the transport speed) denoted $f(t, x; \mathcal{X})$ is regulated by the presence of the axons and the cancer cells. It takes the following form

$$f(t, x; \mathcal{X}) := \pi(x) [1 - \beta(x)\rho(A_1(t)) + \delta(x)A_2(t)] + \eta \left( x, \frac{N_c(t)}{N(t)} \right) \quad (1)$$

where

* $\pi(x)$ represents a basal amplitude for the speed of disease progression. Since we expect the transformation of healthy cells to be slow at early stages, we consider the transfer to be almost negligible for $x \ll 0$. This transport is then expected to increase up to a plateau observed during the early PDAC stage. Moreover, we naturally assume that at the boundary of the phenotype axis $(+\infty$ or $L)$, the function $\pi$ vanishes.

* The basal amplitude is modulated by the presence of the sympathetic axons that slow down the disease progression (cf. Guillot et al., 2022). The function $\beta$ modulates the maximum rate of the regulation depending on the phenotype variable $x$. The function $\rho$ modulates the basal amplitude negatively because of the inhibiting mechanism of $A_1$, under the assumption that a sufficient density of sympathetic axons is required to have an impact on the cancer progression process. Thus, the function $\rho$ vanishes for small values of $A_1$. Details about the mathematical expression and the parameters of the function $\rho$ are given in Section 3.1.

* The basal amplitude is also regulated by the presence of sensory axons that speed up the disease progression (cf. Saloman et al., 2016). The function $\delta$ modulates the maximum rate of the regulation depending on the phenotype variable $x$.

* Moreover, we assume that cancer cells in the pancreas are inducing an acceleration of the cancer progression of the healthy cells in the late stages of the disease. This assumption is formalizing the crosstalk between PDAC cells and other cells where the tumorigenic environment facilitates the ability of cancer cells to survive and proliferate to the detriment of healthy cells (cf. Alexander and Cukierman, 2016). In order to model this response, we introduce the function $\eta$ which contributes to the cancer progression acceleration. This function takes as variables the ratio $N_c(t)/N(t)$ which quantifies the presence and the proportion of pathological cells in the system and the phenotype variable $x$. Also, we assume that $\eta$ is compactly supported and its support is located in a neighborhood corresponding to the healthy cells on the phenotype axis.

An additional assumption in our model is that the cancer progression is non-reversible. For instance, once cells start to progress towards a cancerous state, they cannot recover, meaning that they cannot have an healthy phenotype state later in time. This implies that the transport term $f$ remains non-negative.

**Growth term.** In the model, the growth of the cancer cell densities is modeled by a logistic-type growth which is given by the growth term $g$:

$$g(t, x; \mathcal{X}) := r(x)Q(t, x) \left( 1 - \frac{N(t)}{\tau_C} - \mu_1 A_1(t) + \mu_2 A_2(t) \right). \quad (2)$$

The function $r$ of the phenotype variable $x$ is the basal growth rate of the proliferating cells. Since growth starts at pre-cancerous stages and accelerate around cancerous stages, (cf. Klein et al., 2002), the support of the function $r$ is located in its neighborhood on the phenotype axis. This growth process is regulated by the axons (cf. Biankin et al., 2012). On the one hand, the presence of the sympathetic axons inhibits or promotes the growth of cancer cells depending on the sign of $\mu_1$ which is the amplitude of the modulation (cf. Renz et al., 2018a and Guillot et al., 2022) This mechanism is modeled by the term $-\mu_1 A_1(t)$, if $\mu_1 > 0$, then the sympathetic axons play an inhibitory role, whereas if $\mu_1 < 0$, the sympathetic axons play a promoting role on tumor growth. On the other hand, the presence of sensory axons amplifies the growth of cancer cells (cf. Banh et al., 2020), and this is modeled by the term $+\mu_2 A_2(t)$ where $\mu_2 \ge 0$. Moreover, the saturation rate of the cells densities is linked to the parameters $\tau_C$ (the carrying capacities), $\mu_1$ and $\mu_2$.

**Sympathetic axons growth dynamics.** One interesting dynamics for the sympathetic axons is that a small increase of the density of these axons is observed early in precancerous stage of cancer progression process, whereas once cancer is established the density of sympathetic
* $\pi(x)$ represents a basal amplitude for the speed of disease progression. Since we expect the transformation of healthy cells to be slow at early stages, we consider the transfer to be almost negligible for $x \ll 0$. This transport is then expected to increase up to a plateau observed during the early PDAC stage. Moreover, we naturally assume that at the boundary of the phenotype axis $(+\infty$ or $L)$, the function $\pi$ vanishes.

* The basal amplitude is modulated by the presence of the sympathetic axons that slow down the disease progression (cf. [Guillot et al., 2022]). The function $\beta$ modulates the maximum rate of the regulation depending on the phenotype variable $x$. The function $\rho$ modulates the basal amplitude negatively because of the inhibiting mechanism of $A_1$, under the assumption that a sufficient density of sympathetic axons is required to have an impact on the cancer progression process. Thus, the function $\rho$ vanishes for small values of $A_1$.

* The basal amplitude is also regulated by the presence of sensory axons that speed up the disease progression (cf. [Saloman et al., 2016]). The function $\delta$ modulates the maximum rate of the regulation depending on the phenotype variable $x$.

* Moreover, we assume that cancer cells in the pancreas are inducing an acceleration of the cancer progression of the healthy cells in the late stages of the disease. This assumption is formalizing the crosstalk between PDAC cells and other cells where the tumorigenic environment facilitates the ability of cancer cells to survive and proliferate to the detriment of healthy cells (cf. [Alexander and Cukierman, 2016]). In order to model this response, we introduce the function $\eta$ which contributes to the cancer progression acceleration. This function takes as variables the ratio $N_c(t)/N(t)$ which quantifies the presence and the proportion of pathological cells in the system and the phenotype variable $x$. Also, we assume that $\eta$ is compactly supported and its support is located in a neighbourhood corresponding to the healthy cells on the phenotype axis.

An additional assumption in our model is that the cancer progression is non-reversible. For instance, once cells start to progress towards a cancerous state, they cannot recover, meaning that they cannot have an healthy phenotype state later in time. This implies that the transport term $f$ remains non-negative.

**Growth term.** In the model, the growth of the cancer cell densities is modelled by a logistic-type growth which is given by the growth term $g$:

$$g(t, x; \mathcal{X}) := r(x)Q(x) \left( 1 - \frac{N(t)}{\tau_C} - \mu_1 A_1(t) + \mu_2 A_2(t) \right). \quad (2)$$

The function $r$ of the phenotype variable $x$ is the basal growth rate of the proliferating cells. Since growth starts at pre-cancerous stages and accelerate around cancerous stages, (cf. [Klein et al., 2002]), the support of the function $r$ is located in its neighbourhood on the phenotype axis. This growth process is regulated by the axons (cf. [Biankin et al., 2012]). On the one hand, the presence of the sympathetic axons inhibits or promotes the growth of cancer cells depending on the sign of $\mu_1$ which is the amplitude of the modulation (cf. [Renz et al., 2018a] and [Guillot et al., 2022]) This mechanism is modelled by the term $-\mu_1 A_1(t)$, if $\mu_1 > 0$, then the sympathetic axons play an inhibitory role, whereas if $\mu_1 < 0$, the sympathetic axons play a promoting role on tumor growth. On the other hand, the presence of sensory axons amplifies the growth of cancer cells (cf. [Banh et al., 2020]), and this is modelled by the term $+\mu_2 A_2(t)$ where $\mu_2 \ge 0$. Moreover, the saturation rate of the cells densities is linked to the parameters $\tau_C$ (the carrying capacities), $\mu_1$ and $\mu_2$.

**Sympathetic axons growth dynamics.** One interesting dynamics for the sympathetic axons is that a small increase of the density of these axons is observed early in precancerous stage of cancer progression process, whereas once cancer is established the density of sympathetic axons is reduced (see the biological data in Figure 9 in Section 4.1). Hence, in order to model the time evolution of the sympathetic axons, we use a logistic law with an Allee effect. We denote by $\theta$ the function that enables the Allee effect. The function $\theta$ takes as argument the ratio $N(t)/N(0)$ and it allows to change the dynamics of the sympathetic axons evolution during disease progression. Starting from the healthy state, $N(t) = N(0)$ at least from small time $t$ i.e. the total concentration of cells remains constant, then taking $\theta(1) < A_1(0)$ and the amount of sympathetic axons increases. Moreover, $\theta$ is an increasing function such that if $N(t) > N(0)$ i.e. there are proliferative cells in the system, then at some point $\theta(N(t)/N(0)) > 1$ and the amount of sympathetic axons decreases. Thus the dynamics of the sympathetic axons are given by the following differential equation:

$$\frac{d}{dt}A_1(t) = r_{A_1}A_1(t) \left( \frac{A_1(t)}{\theta \left( \frac{N(t)}{N(0)} \right)} - 1 \right) (1 - A_1(t)), \quad \text{Logistic law with Allee effect}$$

where $r_{A_1} > 0$ is the growth rate of sympathetic axons.

**Sensory axons growth dynamics.** The remodeling of the sensory axons start at cancerous stages (see the biological data in Figure 10 in Section 4.1). Hence, one natural way to model the time evolution of sensory axons is to use a logistic law with a growth rate which is modulated by the presence of cells with non-negative phenotype values.
axons is reduced (see the biological data in Fig. 9 in Section 4.1). Hence, in order to model the time evolution of the sympathetic axons, we use a logistic law with an Allee effect. We denote by $\theta$ the function that enables the Allee effect. The function $\theta$ takes as argument the ratio $N(t)/N(0)$ and it allows to change the dynamics of the sympathetic axons evolution during disease progression. Starting from the healthy state, $N(t) = N(0)$ at least from small time $t$ i.e. the total concentration of cells remains constant, then taking $\theta(1) < A_1(0)$ and the amount of sympathetic axons increases. Moreover, $\theta$ is an increasing function such that if $N(t) > N(0)$ i.e. there are proliferative cells in the system, then at some point $\theta(N(t)/N(0)) > 1$ and the amount of sympathetic axons decreases. Details of the properties of the function $\theta$ are given in Hypothesis 2.4 and its mathematical expression is given in Section 3.1. Thus the dynamics of the sympathetic axons are given by the following differential equation:

$$ \frac{d}{dt} A_1(t) = r_{A_1} A_1(t) \underbrace{\left\{ \frac{A_1(t)}{\theta \left( \frac{N(t)}{N(0)} \right)} - 1 \right\}}_{\text{Logistic law with Allee effect}} (1 - A_1(t)), $$

where $r_{A_1} > 0$ is the growth rate of sympathetic axons.

**Sensory axons growth dynamics.** The remodeling of the sensory axons start at cancerous stages (see the biological data in Fig. 10 in Section 4.1). Hence, one natural way to model the time evolution of sensory axons is to use a logistic law with a growth rate which is modulated by the presence of cells with non-negative phenotype values. We denote $r_{A_2}$ the increasing function taking as argument $N_c(t)/N(t)$ which model the growth rate of sensory axons. The dynamics of the sensory axons are given by the following differential equation:

$$ \frac{d}{dt} A_2(t) = \underbrace{r_{A_2} \left( \frac{N_c(t)}{N(t)} \right)}_{\text{Impact of cancer cells on growth}} \underbrace{A_2(t)(1 - A_2(t))}_{\text{Logistic law}}. $$

For instance, if there is no proliferative cells, then

$$ r_{A_2} \left( \frac{N_c(t)}{N(t)} = 0 \right) = 0 $$

and the amount of sensory axons remains constant. As soon as $N_c(t)/N(t) > 0$ which implies $r_{A_2}(N_c(t)/N(t)) > 0$ then the sensory axons follow the logistic law. Since we assume that $r_{A_2}$ is monotonous, then the more cancer cells there are in the system, the faster the growth of sensory axon is. Details of the properties of the function $r_{A_2}$ are given in Hypothesis 2.4 and its mathematical expression is given in Section 3.1.

**Complete dynamical system.** The system that mathematically formalizes the impact of axons on tumor progression consists in a partial differential equation (PDE) for cell dynamics and two differential equations for axon dynamics. It couples a growth-transport equation for $Q$ with two ODEs with non-local terms for $A_1$ and $A_2$:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial_t Q(t, x) + \partial_x [f(t, x; \mathcal{X}) Q(t, x)]}{\partial t} = g(t, x; \mathcal{X}), \quad t > 0, \quad x \in \Omega \\[1.5em] \displaystyle \frac{d}{dt} A_1(t) = r_{A_1} A_1(t) \left( \frac{A_1(t)}{\theta \left( \frac{N(t)}{N(0)} \right)} - 1 \right) (1 - A_1(t)), \quad t > 0, \\[1.5em] \displaystyle \frac{d}{dt} A_2(t) = r_{A_2} \left( \frac{N_c(t)}{N(t)} \right) A_2(t) (1 - A_2(t)), \quad t > 0. \end{array} \right. \quad (3) $$

The system (3) is completed with the following initial conditions:

$$ Q(0, x) = Q_0(x) \text{ for } x \in \Omega, \quad A_1(0) \cong A_1^{eq} \quad \text{and} \quad A_2(0) = A_2^{eq}, $$

where $A_1^{eq} > 0$ corresponds to the average density of sympathetic axons in a healthy pancreas which has been normalized and $0 < A_1^{eq} \ll 1$ corresponds to the average density a sensory axons in a healthy pancreas also normalized. For modeling purposes, if there is no observations showing the presence of sensory axons in a healthy pancreas, we still consider that there is a negligible amount of sensory axons in the system at the initial state. Otherwise, the growth dynamics could not take place since $A_2(t) = 0$ is an unstable steady-state. One point to note is that if we

consider initially that the pancreas is essentially composed of healthy cells then one way to cope with this assumptions is to assume that $Q_0$ is compactly supported and that the following holds on the support of the initial datum $\text{supp}(Q_0) \subset \Omega/\mathbb{R}_+$. Finally, we add the following boundary condition on the PDE

$$ f(t, x; \mathcal{X}) Q(t, x) = 0, \quad x \in \partial \Omega, $$

which ensures Dirichlet boundary conditions and allows us to neglect any processes acting outside of our domain. This boundary condition is a strong assumption which is sufficient to prove the well-posedness of the model and the fact that no mass is lost on the boundary of the domain. It is ensured by the compact support of the transport term $f$ (cf. Hypothesis 2.2).

## 2.2. Well-posedness of the model

We denote by $\mathcal{X} := (Q, A_1, A_2)$ the solution of the system (3). We note that $\mathcal{X}$ is solution of a non-conservative system that is a particular case of the following system:

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial_t Q(t, x) + f(t, x; \mathcal{X}) \partial_x Q(t, x) + c(t, x; \mathcal{X}) Q(t, x)}{\partial t} = 0 & \text{for } t > 0, x \in \Omega, \quad (a) \\[1.5em] \displaystyle f(t, x; \mathcal{X}) Q(t, x) = 0 & \text{for } x \in \partial \Omega, t \ge 0, \quad (b) \\[1.5em] \displaystyle Q(0, x) = Q_0(x) & \text{for } x \in \Omega, \quad (c) \\[1.5em] \displaystyle \left\{ \begin{array}{l} \frac{d}{dt} A_1(t) = r_{A_1} A_1(t) \left( \frac{A_1(t)}{\theta \left( \frac{N(t)}{N(0)} \right)} - 1 \right) (1 - A_1(t)) \\ A_1(0) = A_1^0 \in (0, 1), \end{array} \right. & \text{for } t > 0, \quad (d) \\[1.5em] \displaystyle A_2(t) = r_{A_2} \left( \frac{N_c(t)}{N(t)} \right) A_2(t) (1 - A_2(t)) & \text{for } t > 0, \quad (e) \\[1.5em] \displaystyle A_2(0) = A_2^0 \in (0, 1), & \text{for } t > 0, \quad (f) \\[1.5em] \displaystyle \frac{d}{dt} A_2(t) = r_{A_2} \left( \frac{N_c(t)}{N(t)} \right) A_2(t) (1 - A_2(t)), & \text{for } t > 0, \quad (g) \end{array} \right. \quad (4) $$

where $N(t) := \int_{\Omega} Q(t, x)dx$ and $N_c(t) := \int_{\Omega} \psi(x)Q(t, x)dx$ such that $\psi$ is given and nonnegative and $\text{supp}(\psi) \subset \Omega \cap \mathbb{R}_+$ and that $\|\psi Q(t)\|_{L^1(\Omega)} \le \|Q(t)\|_{L^1(\Omega)}$ for $t \ge 0$. We also introduce the function $g$ which allows to rewrite Eq. (4) in a conservative form:

$$ g(t, x; \mathcal{X}) := c(t, x; \mathcal{X}) - \partial_x f(t, x; \mathcal{X}), \quad (5) $$

implying that

$$ \partial_t Q(t, x) + \partial_x (f(t, x; \mathcal{X}) Q(t, x)) + g(t, x; \mathcal{X}) Q(t, x) = 0 \quad \text{for } t > 0, \quad x \in \Omega. $$

In the following, we study the well-posedness of system (4) and start with the different hypotheses required for that.

**Hypothesis 2.1 (Initial Condition).** Assume $Q_0 \in C^1(\Omega)$ is nonnegative such that

$$ 0 < C(N(0)) < \int_{\Omega} Q_0(x) dx = \|Q_0\|_{L^1(\Omega)} < C(\tau_c) < \infty $$

with $C(N(0))$ and $C(\tau_c)$ positive constants.

Hence, we denote $\mathcal{P}$ the following set:

$$ \mathcal{P}(T) := \left\{ (Q, A_1, A_2) \in C^1([0, T]; C^1(\Omega) \cap \mathcal{L}^1(\Omega)) \times C^1([0, T]) \times C^1([0, T]) \text{ such that} \right. \\ \left. C(N(0)) \le \int_{\Omega} Q(t, x) dx \le C(\tau_c) \text{ and } 0 \le A_1(t) \le 1, \ i = 1, 2, \ t \in [0, T] \right\}. \quad (6) $$

We also assume the following conditions on the functions given by the transport term and the growth term. Given $\mathcal{X} \in \mathcal{P}(T)$, for $T > 0$, we have the following hypotheses:

**Hypothesis 2.2 (Transport Term).** Assume that $\mathcal{X}_i \in \mathcal{P}(T)$, for $T > 0$ for $i = 1, 2$ then

$$ f(\cdot; \mathcal{X}) \in C([0, T]; C_c^1(\Omega)) \text{ with } |\partial_x f(t, x; \mathcal{X})| < \infty, \ f(t, x; \mathcal{X}) \ge 0 \quad \text{for } t \in [0, T] \text{ and } x \in \Omega, $$
We denote $r_{A_2}$ the increasing function taking as argument $N_c(t)/N(t)$ which model the growth rate of sensory axons. The dynamics of the sensory axons are given by the following differential equation:

$$ \frac{d}{dt} A_2(t) = \underbrace{r_{A_2}\left(\frac{N_c(t)}{N(t)}\right)}_{\text{Impact of cancer cells on growth}} \underbrace{A_2(t)(1 - A_2(t))}_{\text{Logistic law}}. $$

For instance, if there is no proliferative cells, then

$$ r_{A_2}\left(\frac{N_c(t)}{N(t)}\right) = 0 $$

and the amount of sensory axons remains constant. As soon as $N_c(t)/N(t) > 0$ which implies $r_{A_2}(N_c(t)/N(t)) > 0$ then the sensory axons follow the logistic law. Since we assume that $r_{A_2}$ is monotonous, then the more cancer cells there are in the system, the faster the growth of sensory axon is.

**Complete dynamical system.** The system that mathematically formalizes the impact of axons on tumor progression consists in a partial differential equation (PDE) for cell dynamics and two differential equations for axon dynamics. It couples a growth-transport equation for $Q$ with two ODEs with non-local terms for $A_1$ and $A_2$:

$$ \left\{ \begin{array}{l} \partial_t Q(t,x) + \partial_x [f(t,x; \mathcal{X}) Q(t,x)] = g(t,x; \mathcal{X}), \quad t > 0, x \in \Omega \\[1em] \frac{d}{dt} A_1(t) = r_{A_1} A_1(t) \left( \frac{A_1(t)}{\theta\left(\frac{N(t)}{N(0)}\right)} - 1 \right) (1 - A_1(t)), \quad t > 0, \\[1em] \frac{d}{dt} A_2(t) = r_{A_2} \left( \frac{N_c(t)}{N(t)} \right) A_2(t) (1 - A_2(t)), \quad t > 0. \end{array} \right. \quad (3) $$

The system (3) is completed with the following initial conditions:

$$ Q(0, x) = Q_0(x) \text{ for } x \in \Omega, \quad A_1(0) \cong A_1^{eq} \quad \text{and} \quad A_2(0) = A_2^0, $$

where $A_1^{eq} > 0$ corresponds to the average density of sympathetic axons in a healthy pancreas which has been normalized and $0 < A_2^0 \ll 1$ corresponds to the average density a sensory axons in a healthy pancreas also normalized. For modeling purposes, if there is no observations showing the presence of sensory axons in a healthy pancreas, we still consider that there is a negligible amount of sensory axons in the system at the initial state. Otherwise, the growth dynamics could not take place since $A_2(t) = 0$ is an unstable steady-state. One point to note is that if we consider initially that the pancreas is essentially composed of healthy cells then one way to cope with this assumptions is to assume that $Q_0$ is compactly supported and that the following holds on the support of the initial datum $\text{supp}(Q_0) \subset \Omega/\mathbb{R}_+$. Finally, we add the following boundary condition on the PDE

$$ f(t, x; \mathcal{X}) Q(t, x) = 0, x \in \partial\Omega, $$

which ensures Dirichlet boundary conditions and allows us to neglect any processes acting outside of our domain. This boundary condition is a strong assumption which is sufficient to prove the well-posedness of the model and the fact that no mass is lost on the boundary of the domain. It is ensured by the compact support of the transport term $f$ (cf. Hypothesis 2.2).

## 2.2 Well-posedness of the model

We denote by $\mathcal{X} := (Q, A_1, A_2)$ the solution of the system (3). We note that $\mathcal{X}$ is solution of a non-conservative system that is a particular case of the following system :

$$ \partial_t Q(t, x) + f(t, x; \mathcal{X}) \partial_x Q(t, x) + c(t, x; \mathcal{X}) Q(t, x) = 0 \quad \text{for } t > 0, x \in \Omega, \quad (4a) $$

$$ f(t, x; \mathcal{X}) Q(t, x) = 0 \quad \text{for } x \in \partial\Omega, t \ge 0, \quad (4b) $$

$$ Q(0, x) = Q_0(x) \quad \text{for } x \in \Omega, \quad (4c) $$

$$ \frac{d}{dt} A_1(t) = r_{A_1} A_1(t) \left( \frac{A_1(t)}{\theta\left(\frac{N(t)}{N(0)}\right)} - 1 \right) (1 - A_1(t)) \quad \text{for } t > 0, \quad (4d) $$

$$ A_1(0) = A_1^0 \in (0, 1), \quad (4e) $$

$$ \frac{d}{dt} A_2(t) = r_{A_2} \left( \frac{N_c(t)}{N(t)} \right) A_2(t) (1 - A_2(t)) \quad \text{for } t > 0, \quad (4f) $$

$$ A_2(0) = A_2^0 \in (0, 1), \quad (4g) $$
$$f(t, x; \mathcal{X}) = \partial_x f(t, x; \mathcal{X}) = 0 \text{ for } t \in [0, T] \text{ and } x \in \partial\Omega$$

and there exist constants $C_l(f) > 0$, $C_l(\partial_x f) > 0$ such that

$$\|f(t; \mathcal{X}_1) - f(t; \mathcal{X}_2)\|_{L^\infty(\Omega)} \le C_l(f)\|\mathcal{X}_1 - \mathcal{X}_2\|_p \quad \text{and} \\ \|\partial_x(f(t; \mathcal{X}_1) - f(t; \mathcal{X}_2))\|_{L^\infty(\Omega)} \le C_l(\partial_x f)\|\mathcal{X}_1 - \mathcal{X}_2\|_p.$$

**Hypothesis 2.3 (Growth Term).** Assume that $\mathcal{X}_l \in P(T)$, for $T > 0$ for $i = 1, 2$ then

$$g(\cdot; \mathcal{X}) \in C([0, T]; C_c^1(\Omega)) \text{ with } |g(t, x; \mathcal{X})| < \infty, |\partial_x g(t, x; \mathcal{X})| < \infty \\ \text{for } t \in [0, T] \text{ and } x \in \Omega$$

and there exists a constant $C_l(g) > 0$ such that

$$\|g(t; \mathcal{X}_1) - g(t; \mathcal{X}_2)\|_{L^\infty(\Omega)} \le C_l(g)\|\mathcal{X}_1 - \mathcal{X}_2\|_p.$$

Moreover, assume that the growth term is a logistic-type growth term, i.e

$$-g(t, x; \mathcal{X}) := r(x) \left( h \left( A_1(t), A_2(t) \right) - \int_{\Omega} Q(t, y) dy \right)$$

where there exist constants $0 \le r_- < r_+ < \infty$ and $C(\tau_c) > 0$ such that

$$r_- \le r(x) \le r_+ \quad \text{for } x \in \Omega$$

and there exists a small perturbation $0 < \varepsilon$ such that

$$C(\tau_c) - \varepsilon \le h(x, y) \le C(\tau_c) \quad \text{for } (x, y) \in [0, 1]^2.$$

**Hypothesis 2.4 (Coupling Axons and Cell Densities).** Assume $\theta \in C(\mathbb{R}_+)$ and $r_{A_2} \in C(\mathbb{R}_+)$ such that

$$r_{A_2}(x) \ge 0, \quad \forall x \in \mathbb{R}$$

and that for any compact set $K \subset \mathbb{R}_+$, there exist a constant $C_l(\theta) > 0$ and a constant $C_l(r_{A_2}) > 0$ such that $\forall (x, y) \in K \times K$

$$|\theta(x) - \theta(y)| \le C_l(\theta)|x - y| \quad \text{and} \quad |r_{A_2}(x) - r_{A_2}(y)| \le C_l(r_{A_2})|x - y|.$$

Moreover, the image of $\theta$ satisfies

$$\mathrm{Im}(\theta) \subset [\theta_-, \theta_+] \text{ with } \theta_- \in (0, A_1^0), \theta_+ > 1.$$

**Remark 1.** The **Hypothesis 2.1** allows to consider an initial state of healthy cells for the system. As for **Hypotheses 2.2 and 2.3**, they give sufficient conditions on the tumor progression term (the transport term $f$) and on the proliferation term (the reaction term $c$) in order to ensure the well-posedness. The **Hypothesis 2.4** states some regularity assumptions on the terms formalizing the effect of the cells on the axons densities as well as additional conditions in order to have a coupling term biologically relevant (cf. Section 2.1).

**Remark 2.** In our model,

$$f(t, x; \mathcal{X}) = \pi(x)(1 - \beta(x)\rho(A_1(t)) + \delta(x)A_2(t)) + \eta \left( x, \frac{N_c(t)}{N(t)} \right)$$

and

$$c(t, x; \mathcal{X}) = \partial_x f(t, x; \mathcal{X}) - r(x) \left( 1 - \frac{N(t)}{\tau_c} - \mu_1 A_1(t) + \mu_2 A_2(t) \right).$$

Hence the conditions on $f$ and $c$ stated in **Hypotheses 2.2 and 2.3** are enforced when the following detailed assumptions are fulfilled:

* Let $\pi \in C_c^1(\Omega)$, $\beta \in C^1(\Omega)$, $\delta \in C^1(\Omega)$ and $\eta \left( \cdot, \frac{N_c(t)}{N(t)} \right) \in C_c^1(\Omega)$ such that

$$\sup_{x \in \Omega} |\pi'(x)| + |\beta'(x)| + |\delta'(x)| + |\eta'(x)| \le C$$

for a constant $C > 0$.

* Let $r \in C_c^1(\Omega)$ and Lipschitz continuous.

More details about these functions will be given in Section 3.1.

**Theorem 1.** Assume that the **Hypothesis 2.1, 2.2, 2.3 and 2.4** are satisfied, then the system (4) admits a unique solution

$$(Q, A_1, A_2) \in C^1([0, T]; C^1(\Omega) \cap L^1(\Omega)) \times C^1([0, T]) \times C^1([0, T]) \text{ for all } T > 0$$

such that

$$\forall s \in [0, T] \quad 0 \le A_i(s) \le 1, \quad i = 1, 2,$$

and there exists a positive constant $C(\tau_c) > 0$

$$\forall s \in [0, T] \quad \int_{\Omega} Q_0(x) dx \le \int_{\Omega} Q(s, x) dx \le C(\tau_c).$$

**Proof.** The system (4) couples an transport-reaction partial differential equation (PDE) for the cancer progression with two differential equations for the axons densities. Axons dynamic is governed by non-local terms that depend on the solution of the PDE. Inversely, the dynamic of tumor progression modeled by the transport term and the reaction term is also governed by the axons. The proof of the well-posedness of the solution of this non-linear system relies on the contraction mapping theorem.

Consider $\mathcal{X} = (Q, A_1, A_2) \in C^1([0, T]; C^1(\Omega) \cap L^1(\Omega)) \times C^1([0, T]) \times C^1([0, T])$ given. Define $N(s) := \int_{\Omega} Q(s, x) dx$ and $N_c(s) := \int_{\Omega} \psi(x) Q(s, x) dx$ such that $0 < N(0) \le N(s) \le C(\tau_c)$ for $s \in [0, T]$. We introduce the linear system associated to (4):

$$\left\{
\begin{array}{ll}
\displaystyle \partial_t u(s, x) + f(s, x; \mathcal{X}) \partial_x u(s, x) + c(s, x; \mathcal{X}) u(s, x) = 0, & \text{for } s \in (0, T), x \in \Omega, \\
\\
\displaystyle f(s, x) u(s, x) = 0, & \text{for } x \in \partial\Omega, s \in (0, T), \\
\\
\displaystyle u(0, x) = Q_0(x), & \text{for } x \in \Omega, \\
\\
\begin{array}{ll}
\displaystyle \frac{d}{ds} \widetilde{A}_1(s) = r_{A_1} \widetilde{A}_1(s) \left( \frac{\widetilde{A}_1(s)}{\theta \left( \frac{N(s)}{N(0)} \right)} - 1 \right) (1 - \widetilde{A}_1(s)), & \widetilde{A}_1(0) = A_1^0, \\
\\
\displaystyle \frac{d}{ds} \widetilde{A}_2(s) = r_{A_2} \left( \frac{N_c(s)}{N(s)} \widetilde{A}_2(s) (1 - \widetilde{A}_2(s)), \quad \widetilde{A}_2(0) = A_2^0 \right).
\end{array}
\end{array}
\right.
\quad (7)$$

We denote

$$B := \left\{ \begin{array}{l} Q \in C^1([0, T]; C^1(\Omega) \cap L^1(\Omega)) \\ \le \int_{\Omega} Q(s, x) dx \le C(\tau_c) \text{ for } s \in [0, T] \end{array} \right\}$$

where $C(N(0))$ is a constant such that $0 < C(N(0)) \le N(0)$, the norm on $B$ is given by

$$\|Q\|_B := \sup_{0 \le s \le T} \|Q(s)\|_{L^1(\Omega)}.$$

and $S$ is the following mapping:

$$S: B \times C^1([0, T]) \times C^1([0, T]) \to B \times C^1([0, T]) \times C^1([0, T]), \quad (Q, A_1, A_2) \mapsto (u, \widetilde{A}_1, \widetilde{A}_2). \quad (8)$$

The norm associated to the set $P(T)$ (we use the notation $P = P(T)$ in the following for the sake of clarity) defined in (6) or by

$$P := B \times \left\{ A \in C^1([0, T]) \mid 0 \le A(t) \le 1, \text{ for } t \in [0, T] \right\}^2,$$

is given by

$$\|(Q, A_1, A_2)\|_P := \left( \|Q\|_B + \sup_{0 \le s \le T} |A_1(s)| + \sup_{0 \le s \le T} |A_2(s)| \right).$$

The proof of the well-posedness of the system (4) works as follows. First, we prove that the linear system (7) associated to (4) admits a unique solution. Second, we prove that the map $S$ is a contraction on the set $P(T)$ for $T$ small enough. Third, we use a bootstrap argument to conclude the existence of a solution of (4). The final step is the uniqueness result. All these steps relies on Lemmas detailed in Appendix A.

**Step 1.** For $\mathcal{X}' = (Q, A_1, A_2) \in P$ given, thanks to Lemmas 1, 2 and 3, the linear problem (7) admits a unique solution
where $N(t) := \int_{\Omega} Q(t, x)dx$ and $N_c(t) := \int_{\Omega} \psi(x)Q(t, x)dx$ such that $\psi$ is given and nonnegative and $\text{supp}(\psi) \subset \Omega \cap \mathbb{R}_+$ and that $\|\psi Q(t)\|_{L^1(\Omega)} \le \|Q(t)\|_{L^1(\Omega)}$ for $t \ge 0$. We also introduce the function $g$ which allows to rewrite the equation (4a) in a conservative form:

$$g(t, x; \mathcal{X}) := c(t, x; \mathcal{X}) - \partial_x f(t, x; \mathcal{X}), \qquad (5)$$

implying that

$$\partial_t Q(t, x) + \partial_x (f(t, x; \mathcal{X}) Q(t, x)) + g(t, x; \mathcal{X}) Q(t, x) = 0 \quad \text{for } t > 0, x \in \Omega.$$

In the following, we study the well-posedness of system (4) and start with the different hypotheses required for that.

**Hypothesis 2.1 (Initial condition).** Assume $Q_0 \in C^1(\Omega)$ is nonnegative such that

$$0 < C(N(0)) < \int_{\Omega} Q_0(x)dx = \|Q_0\|_{L^1(\Omega)} < C(\tau_c) < \infty$$

with $C(N(0))$ and $C(\tau_c)$ positive constants.

Hence, we denote $\mathcal{P}$ the following set:

$$\mathcal{P}(T) := \left\{ (Q, A_1, A_2) \in C^1([0, T]; C^1(\Omega) \cap L^1(\Omega)) \times C^1([0, T]) \times C^1([0, T]) \text{ such that} \right. \\ \left. C(N(0)) \le \int_{\Omega} Q(t, x)dx \le C(\tau_c) \text{ and } 0 \le A_i(t) \le 1, i = 1, 2, t \in [0, T] \right\}. \qquad (6)$$

We also assume the following conditions on the functions given by the transport term and the growth term. Given $\mathcal{X} \in \mathcal{P}(T)$, for $T > 0$, we have the following hypotheses:

**Hypothesis 2.2 (Transport term).** Assume that $\mathcal{X}_i \in \mathcal{P}(T)$, for $T > 0$ for $i = 1, 2$ then

$$f(\cdot; \mathcal{X}) \in C([0, T]; C_c^1(\Omega)) \text{ with } |\partial_x f(t, x; \mathcal{X})| < \infty, f(t, x; \mathcal{X}) \ge 0$$

for $t \in [0, T]$ and $x \in \Omega$,

$$f(t, x; \mathcal{X}) = \partial_x f(t, x; \mathcal{X}) = 0 \text{ for } t \in [0, T] \text{ and } x \in \partial\Omega$$

and there exist constants $C_l(f) > 0$, $C_l(\partial_x f) > 0$ such that

$$\|f(t; \mathcal{X}_1) - f(t; \mathcal{X}_2)\|_{L^\infty(\Omega)} \le C_l(f) \|\mathcal{X}_1 - \mathcal{X}_2\|_{\mathcal{P}} \quad \text{and} \quad \|\partial_x(f(t; \mathcal{X}_1) - f(t; \mathcal{X}_2))\|_{L^\infty(\Omega)} \le C_l(\partial_x f) \|\mathcal{X}_1 - \mathcal{X}_2\|_{\mathcal{P}}.$$

**Hypothesis 2.3 (Growth term).** Assume that $\mathcal{X}_i \in \mathcal{P}(T)$, for $T > 0$ for $i = 1, 2$ then

$$g(\cdot; \mathcal{X}) \in C([0, T]; C_c^1(\Omega)) \text{ with } |g(t, x; \mathcal{X})| < \infty, |\partial_x g(t, x; \mathcal{X})| < \infty \text{ for } t \in [0, T] \text{ and } x \in \Omega$$

and there exists a constant $C_l(g) > 0$ such that

$$\|g(t; \mathcal{X}_1) - g(t; \mathcal{X}_2)\|_{L^\infty(\Omega)} \le C_l(g) \|\mathcal{X}_1 - \mathcal{X}_2\|_{\mathcal{P}}.$$

Moreover, assume that the growth term is a logistic-type growth term, i.e

$$-g(t, x; \mathcal{X}) := r(x) \left( h(A_1(t), A_2(t)) - \int_{\Omega} Q(t, y)dy \right)$$

where there exist constants $0 \le r_- < r_+ < \infty$ and $C(\tau_c) > 0$ such that

$$r_- \le r(x) \le r_+ \quad \text{for } x \in \Omega$$

and there exists a small perturbation $0 < \varepsilon$ such that

$$C(\tau_c) - \varepsilon \le h(x, y) \le C(\tau_c) \quad \text{for } (x, y) \in [0, 1]^2.$$

**Hypothesis 2.4 (Coupling axons and cell densities).** Assume $\theta \in C(\mathbb{R}_+)$ and $r_{A_2} \in C(\mathbb{R}_+)$ such that

$$r_{A_2}(x) \ge 0, \quad \forall x \in \mathbb{R}$$

and that for any compact set $K \subset \mathbb{R}_+$, there exist a constant $C_l(\theta) > 0$ and a constant $C_l(r_{A_2}) > 0$ such that $\forall (x, y) \in K \times K$

$$|\theta(x) - \theta(y)| \le C_l(\theta)|x - y| \quad \text{and} \quad |r_{A_2}(x) - r_{A_2}(y)| \le C_l(r_{A_2})|x - y|.$$

Moreover, the image of $\theta$ satisfies

$$\mathrm{Im}(\theta) \subset [\theta_-, \theta_+] \text{ with } \theta_- \in (0, A_1^0), \theta_+ > 1.$$
$$ (u, \tilde{A}_1, \tilde{A}_2) \in C^1([0, T]; C^1(\Omega) \cap L^1(\Omega)) \times C^1([0, T]) \times C^1([0, T]). $$

Also, thanks to Lemma 3, choosing $T$ such that

$$ T < \frac{\log \left( \frac{C(\tau_C)}{N(0)} \right)}{\|r\|_{\infty} |C(\tau_C) - N(0)|} \times \frac{C(\tau_C) - N(0)}{C(\tau_C) - C(N(0))} = \frac{\lambda_1}{\|r\|_{\infty} \xi_1} \quad (9) $$

where $\xi_1 \in (N(0); C(\tau_C))$ and $\lambda_1 = \frac{C(\tau_C) - N(0)}{C(\tau_C) - C(N(0))}$ and that

$$ T < \frac{\log \left( \frac{N(0)}{C(N(0))} \right)}{\|r\|_{\infty} |N(0) - C(N(0))|} \times \frac{N(0) - C(N(0))}{C(\tau_C) - C(N(0))} = \frac{\lambda_2}{\|r\|_{\infty} \xi_2} \quad (10) $$

where $\xi_2 \in (C(N(0)); N(0))$ and $\lambda_2 = \frac{N(0) - C(N(0))}{C(\tau_C) - C(N(0))}$ then the solution $u$ of the linear system (7) is such that $u \in B$.

Step 2. In order to prove that $S$ has a fixed point, we introduce the following maps :

$$ \mathcal{F}: \mathcal{P} \to \mathcal{B}, \quad \mathcal{F}(\mathcal{X}) = u, $$

$$ \Gamma_1 : \mathcal{P} \to C^1([0, T]), \quad \Gamma_1(\mathcal{X}) = \tilde{A}_1, $$

$$ \Gamma_2 : \mathcal{P} \to C^1([0, T]), \quad \Gamma_2(\mathcal{X}) = \tilde{A}_2. $$

The contraction property of $S$ is implied by some stability properties on the mappings $\mathcal{F}$, $\Gamma_1$ and $\Gamma_2$. In the following, we denote $\mathcal{X}_i = (Q_i, A_1^i, A_2^i)$ for $i = 1, 2$.

First, using the results in Lemma 1 and the assumptions of the function $\theta$ stated in Hypothesis 2.4, we note that

$$ \sup_{0 \le t \le T} |\Gamma_1(\mathcal{X}_1) - \Gamma_1(\mathcal{X}_2)| \le C \left( r_{A_1}, \frac{1}{(\theta-\bar{\theta})^2} \right) Te^{C(\theta)T} \sup_{0 \le t \le T} \left| \theta \left( \frac{N(t; \mathcal{X}_1)}{N(0)} \right) - \theta \left( \frac{N(t; \mathcal{X}_2)}{N(0)} \right) \right| \\ \le C \left( r_{A_1}, \frac{1}{(\theta-\bar{\theta})^2}, C_i(\theta) \right) Te^{C(\theta)T} \sup_{0 \le t \le T} \|Q_1(t) - Q_2(t)\|_{L^1(\Omega)}, \quad (11) $$

where $0 < C \left( r_{A_1}, \frac{1}{(\theta-\bar{\theta})^2}, C_i(\theta) \right) < \infty$ and $0 < C(\theta) < \infty$ are two constants.

Similarly, using the results in Lemma 2 and the assumptions on the function $r_{A_2}$ stated in Hypothesis 2.4, we note that

$$ \sup_{0 \le t \le T} |\Gamma_2(\mathcal{X}_2) - \Gamma_2(\mathcal{X}_1)| \le C(A_2^0)T \sup_{0 \le t \le T} |r_{A_2} \left( \frac{N_c(t; \mathcal{X}_1)}{N(t; \mathcal{X}_1)} \right) - r_{A_2} \left( \frac{N_c(t; \mathcal{X}_2)}{N(t; \mathcal{X}_2)} \right)|, \\ \le C(A_2^0, C_i(r_{A_2}))T \sup_{0 \le t \le T} \left| \frac{N_c(t; \mathcal{X}_1)N(t; \mathcal{X}_2) - N_c(t; \mathcal{X}_2)N(t; \mathcal{X}_1)}{N(t; \mathcal{X}_1)N(t; \mathcal{X}_2)} \right|, \\ \le C(A_2^0, C_i(r_{A_2}), \frac{1}{C(N(0))^2}, C(\tau_C))T \sup_{0 \le t \le T} \|Q_1(t) - Q_2(t)\|_{L^1(\Omega)}, \quad (12) $$

where $0 < C(A_2^0, C_i(r_{A_2}), \frac{1}{C(N(0))^2}, C(\tau_C)) < \infty$ is a constant. Moreover, one can check that the constant $C(A_2^0, C_i(r_{A_2}), \frac{1}{C(N(0))^2}, C(\tau_C))$ is bounded by above and below by strictly positive constants independent of the initial data $A_2^0$ and $Q_0$ since there exist constants $m_1 > 0$ and $m_2 > 0$ such that

$$ m_1 < A_2^0 < 1 \quad \text{and} \quad m_2 < C(N(0)) < C(\tau_C). $$

Also, thanks to the stability result of Lemma 3, we note that

$$ \sup_{0 \le t \le T} \|\mathcal{F}(\mathcal{X}_1) - \mathcal{F}(\mathcal{X}_2)\|_{L^1(\Omega)} \le C(T)T \|\mathcal{X}_1 - \mathcal{X}_2\|_p. \quad (13) $$

Finally, using the inequalities (11), (12) and (13), we get that there exist two constants $0 < C_1 < \infty$ and $0 < C_2 < \infty$ such that

$$ \|S(\mathcal{X}_1) - S(\mathcal{X}_2)\|_p \le C_1(T)Te^{C_2T} \|\mathcal{X}_1 - \mathcal{X}_2\|_p $$

where, using the notations introduced in Lemmas 1, 2 and 3 and those of Hypotheses 2.1, 2.2, 2.3 and 2.4, we have

• $C_1(T) = C(C_i(r_{A_2}), A_2^0, \frac{1}{C(N(0))^2}, C(\tau_C), r_{A_1}, \frac{1}{(\theta-\bar{\theta})^2}, C_i(\theta), C_i(\partial_x f), \|\partial_x g\|_{\infty}, C_i(g), T)$,

• $C_2 = C(C(\theta), r_+, |C(\tau_C) - C(N(0))|)$.

Moreover, we introduce the new constants $\tilde{C}_1(T)$ and $\tilde{C}_2$ independent of $A_2^0$ and $C(N(0))$ such that

$$ C_1(T) \le \tilde{C}_1 \left( C_i(r_{A_2}), \frac{1}{(m_2)^2}, C(\tau_C), r_{A_1}, \frac{1}{(\theta-\bar{\theta})^2}, C_i(\theta), C_i(\partial_x f), \|\partial_x g\|_{\infty}, C_i(g), T \right) $$

and

$$ C_2 \le \tilde{C}_2 (C(\theta), r_+, |C(\tau_C) - m_2|). $$

Since the map

$$ g: t \in [0, T] \mapsto g(t) := \tilde{C}_1(t)te^{\tilde{C}_2t} $$

is continuous, increasing and $g(0) = 0$, there exists $T_1 > 0$ such that

$$ \tilde{C}_1(T_1)T_1e^{\tilde{C}_2T_1} < 1. \quad (14) $$

Choosing $T_1 > 0$ satisfying (9), (10) and (14) with $T_1 < \bar{T}$ that is

$$ T_1 < \min \left( \frac{\lambda_1}{r_+ C(\tau_C)}, \frac{e^{-\tilde{C}_2 \bar{T}}}{\tilde{C}_1(\bar{T})}, \frac{\lambda_2}{r_+ C(\tau_C)} \right), $$

it follows that $S$ is a strict contraction on $\mathcal{P}$. Consequently, the contraction mapping Theorem implies there exists a unique $\mathcal{X} \in \mathcal{P}$ of (4) on the time interval $[0, T_1]$.

Step 3. Now, we prove that we can extend the solution to $[0, T]$ for any $T > 0$. In order to do so, since $Q(t) \in C^1(\Omega) \cap L^1(\Omega)$ for $0 \le t \le T_1$, we prove that we can repeat the previous arguments to extend the solution to the time interval $[T_1, 2T_1]$. We first note that Lemma 4 ensures that any solution $Q$ of (4) satisfies

$$ N(0) \le \int_{\Omega} Q(t, x) dx \le C(\tau_C). $$

It implies that $N(T_1) \ge N(0)$. Hence, there exists a positive constant $C(N(T_1))$ such that

$$ 0 < C(N(T_1)) < N(T_1) \le C(\tau_C) $$

and that

$$ C(N(T_1)) = N(T_1) - (N(0) - C(N(0))) \frac{C(\tau_C) - N(T_1)}{C(\tau_C) - N(0)}. \quad (15) $$

Moreover, we define the space on which we apply the contraction mapping Theorem as the following

$$ P_1 = B_1 \times \{A \in C^1([T_1, 2T_1]) \mid 0 \le A(t) \le 1, \text{ for } t \in [T_1, 2T_1]\}^2, $$

where

$$ B_1 := \{ Q \in C^1([T_1, 2T_1]; C^1(\Omega) \cap L^1(\Omega)) \mid C(N(T_1)) \\ \le \int_{\Omega} Q(s, x) dx \le C(\tau_C) \text{ for } s \in [T_1, 2T_1] \}. $$

Hence, we can proceed similarly as in Step 1 and Step 2. Thanks to (15), we have

$$ \lambda_1 = \frac{C(\tau_C) - N(0)}{C(\tau_C) - C(N(0))} = \frac{C(\tau_C) - N(T_1)}{C(\tau_C) - C(N(T_1))} $$

and the dependance to the initial condition which appears in (9) does not cause an issue since

$$ \frac{\lambda_1}{r_+ C(\tau_C)} \le \lambda_1 \frac{\log \left( \frac{C(\tau_C)}{N(T_1)} \right)}{\|r\|_{\infty} |C(\tau_C) - N(T_1)|} = \frac{\lambda_1}{\|r\|_{\infty} \xi} $$

where $\xi \in (N(T_1); C(\tau_C))$. As for the dependance to the initial condition which appears in (10), thanks to (15), we have

$$ \lambda_2 = \frac{N(0) - C(N(0))}{C(\tau_C) - C(N(0))} = \frac{N(T_1) - C(N(T_1))}{C(\tau_C) - C(N(T_1))} $$

and

$$ \frac{\lambda_2}{r_+ C(\tau_C)} \le \lambda_2 \frac{\log \left( \frac{N(T_1)}{C(N(T_1))} \right)}{\|r\|_{\infty} |N(T_1) - C(N(T_1))|} = \frac{\lambda_2}{\|r\|_{\infty} \xi} $$

where $\xi \in (C(N(T_1)); N(T_1))$.

Also, thanks to (15) and the fact that $N(0) \le N(T_1)$, the following holds
**Remark 1.** The Hypothesis 2.1 allows to consider an initial state of healthy cells for the system. As for Hypotheses 2.2 and 2.3, they give sufficient conditions on the tumor progression term (the transport term f) and on the proliferation term (the reaction term c) in order to ensure the well-posedness. The Hypothesis 2.4 states some regularity assumptions on the terms formalizing the effect of the cells on the axons densities as well as additional conditions in order to have a coupling term biologically relevant (cf. Section 2.1).

**Remark 2.** In our model,

$$f(t, x; \mathcal{X}) = \pi(x)(1 - \beta(x)\rho(A_1(t)) + \delta(x)A_2(t)) + \eta \left(x, \frac{N_c(t)}{N(t)}\right)$$

and

$$c(t, x; \mathcal{X}) = \partial_x f(t, x; \mathcal{X}) - r(x) \left( 1 - \frac{N(t)}{\tau_c} - \mu_1 A_1(t) + \mu_2 A_2(t) \right).$$

Hence the conditions on f and c stated in Hypotheses 2.2 and 2.3 are enforced when the following detailed assumptions are fulfilled:

* Let $\pi \in C_c^1(\Omega)$, $\beta \in C^1(\Omega)$, $\delta \in C^1(\Omega)$ and $\eta \left(\cdot, \frac{N_c(t)}{N(t)}\right) \in C_c^1(\Omega)$ such that
$$ \sup_{x \in \Omega} |\pi'(x)| + |\beta'(x)| + |\delta'(x)| + |\eta'(x)| \le C $$

for a constant $C > 0$.

* Let $r \in C_c^1(\Omega)$ and Lipschitz continuous.

More details about these functions will be given in Section 3.1.

**Theorem 1.** Assume that the Hypotheses 2.1, 2.2, 2.3 and 2.4 are satisfied, then the system (4) admits a unique solution

$$(Q, A_1, A_2) \in C^1([0, T]; \mathcal{L}^1(\Omega)) \times C^1([0, T]) \times C^1([0, T]), \text{ for all } T > 0$$

such that

$$\forall s \in [0, T] \quad 0 \le A_i(s) \le 1, \quad i = 1, 2,$$

and there exists a positive constant $C(\tau_c) > 0$

$$\forall s \in [0, T] \quad \int_{\Omega} Q_0(x) dx \le \int_{\Omega} Q(s, x) dx \le C(\tau_c).$$

**Sketch of the proof:** The details of the proof of Theorem 1 are postponed in the Appendix A. The system (4) couples an transport-reaction partial differential equation (PDE) for the cancer progression with two differential equations for the axons densities. Axons dynamic is governed by non-local terms that depend on the solution of the PDE. Inversely, the dynamic of tumor progression modelled by the transport term and the reaction term is also governed by the axons. The proof of the well-posedness of the solution of this non-linear system relies on the contraction mapping theorem.

Consider $\mathcal{X} = (Q, A_1, A_2) \in C^1([0, T]; C^1(\Omega) \cap \mathcal{L}^1(\Omega)) \times C^1([0, T]) \times C^1([0, T])$ given. Define $N(s) := \int_{\Omega} Q(s, x) dx$ and $N_c(s) := \int_{\Omega} \psi(x) Q(s, x) dx$ such that $0 < N(0) \le N(s) \le C(\tau_c)$ for $s \in [0, T]$. We introduce the linear system (4):

$$ 
\left\{
\begin{array}{ll}
\displaystyle \partial_s u(s, x) + f(s, x; \mathcal{X}) \partial_x u(s, x) + c(s, x; \mathcal{X}) u(s, x) = 0, & \text{for } s \in (0, T), x \in \Omega, \\
\\
\displaystyle f(s, x) u(s, x) = 0, & \text{for } x \in \partial \Omega, s \in (0, T), \\
\\
\displaystyle u(0, x) = Q_0(x), & \text{for } x \in \Omega, \\
\\
\displaystyle \frac{d}{ds} \widetilde{A}_1(s) = r_{A_1} \widetilde{A}_1(s) \left( \frac{\widetilde{A}_1(s)}{\theta \left( \frac{N(s)}{N(0)} \right)} - 1 \right) (1 - \widetilde{A}_1(s)), & \widetilde{A}_1(0) = A_1^0, \\
\\
\displaystyle \frac{d}{ds} \widetilde{A}_2(s) = r_{A_2} \left( \frac{N_c(s)}{N(s)} \right) \widetilde{A}_2(s) (1 - \widetilde{A}_2(s)), & \widetilde{A}_2(0) = A_2^0.
\end{array}
\right.
\quad (7) $$

We denote

$$B := \left\{ Q \in C^1([0, T]; C^1(\Omega) \cap \mathcal{L}^1(\Omega)) \mid \int_{\Omega} Q_0(x) dx \le \int_{\Omega} Q(s, x) dx \le C(\tau_c) \text{ for } s \in [0, T] \right\}$$

and $S$ the following mapping:

$$S: B \times C^1([0, T]) \times C^1([0, T]) \to B \times C^1([0, T]) \times C^1([0, T]), \quad (Q, A_1, A_2) \mapsto (u, \widetilde{A}_1, \widetilde{A}_2). \quad (8)$$

Hence, we note that the solution of (4) is a fixed point of $S$. The Lemmas 1, 2 and 3 in Appendix A prove the well-posedness of the solution of (7) and give the estimations needed to prove the contraction property of the mapping which implies the existence and uniqueness of the solution of (4).
$$C(N(T_1)) - C(N(0)) = \frac{C(\tau_C) - C(N(0))}{C(\tau_C) - N(0)} (N(T_1) - N(0)) \ge 0.$$

Then, we have that

$$\frac{e^{-\tilde{C}_2 \tilde{T}}}{\tilde{C}_1(\tilde{T})} \le \frac{e^{-\tilde{C}_2 C(N(T_1)) \tilde{T}}}{C_1(\tilde{T}, C(N(T_1)))}.$$

Using Lemmas 1, 2 and 3 and since $2T_1 - T_1 = T_1 < \tilde{T}$, we obtain the strict contraction on the mapping on $P_1$ with the same condition

$$T_1 < \min \left( \frac{\lambda_1}{r_+ C(\tau_C)}, \frac{e^{-\tilde{C}_2 \tilde{T}}}{\tilde{C}_1(\tilde{T})}, \frac{\lambda_2}{r_+ C(\tau_C)} \right).$$

Hence, by iterating, we extend the solution to the full interval $[0, T]$.

**Step 4.** The uniqueness of the solution of the system is direct consequence of the stability properties of Lemmas 1, 2 and 3 and the structural property of our system which implies that $\mathcal{X} = 0$ is a stationary solution. ■

**Remark 3.** The Theorem 1 states the well-posedness on any finite time interval of the solution of (4). However, one can prove the global well-posedness with the same arguments assuming additional regularity assumptions (on the time variable $t \in \mathbb{R}_+$ instead of $t \in [0, T]$) on the functions $\theta$ and $r_{A_2}$ (e.g. no singularities). The main reason allowing the extension of the result about the well-posedness comes from the fact that assuming **Hypothesis 2.1** holds and $0 < A_1^0 < 1$ for $i = 1, 2$ then the trajectories $(\int_{\Omega} Q(\cdot, x)dx, A_1(t), A_2(t))$ remain bounded away from 0. However, finer stability estimates are required to conclude. Since this study focuses on the interactions between cells and axons during the tumorigenesis for the pancreatic cancer, the transient behavior (instead of the asymptotic behavior) is the main point of interest. The result of Theorem 1 is sufficient to pursue this goal.

## 2.3. Qualitative study of the system

The aim of this section is to study how the transport term (1) and the proliferation term (2) of the system (4) impact the final state of the model's dynamics. First, we show that, depending on the assumptions made about the functions $\pi$ (the basal speed of disease progression), $\eta$ (the speed of disease progression due to the presence of cancer cells), $r$ (the basal proliferation rate of cancer cells) and $Q_0$ (the initial distribution of cells with respect to their phenotypes), the different behaviors of the model's dynamics can be summed up in three categories. Secondly, we obtain explicit bounds on the first time of appearance of cells with cancerous phenotypes. This result is obtained under some assumptions on the transport term and relies on the characteristic curves of the equation.

### 2.3.1. Behavior of the model

We recall that our phenotype domain $\Omega \subset \mathbb{R}$ can either be the whole real line or a segment centered in 0. Hence, a distinct separation for the cells with positive and negative phenotypes is assumed. The origin separates the cancer-induced proliferative cells from the healthy one. Moreover, at time $t = 0$ the cells are supposed to be at a non-cancerous and thus at a non-proliferating state. This assumption is formally translated in our model by the following hypothesis on the location on the phenotype axis of the initial distribution $Q_0$ (i.e. its support) and on the support of the proliferation rate $r$:

**Hypothesis 2.5.** Assume $Q_0$ is the initial condition of (4) which fulfills **Hypothesis 2.1**. Moreover, let $Q_0$ be compactly supported and

$$\mathrm{supp}(Q_0) \cap \mathrm{supp}(r) = \emptyset.$$

**Hypothesis 2.5** implies that the population of cells in a healthy pancreas remain constant until the apparition of cancer cells. Moreover, we also assume the following:

### **Hypothesis 2.6.** Let

$$\pi(x) \ge 0, \quad \eta(x) \ge 0, \quad \delta(x) \ge 0 \quad \text{and} \quad 1 - \beta(x)\rho(y) \ge 0 \quad \forall x \in \Omega, \forall y \in [0, 1].$$

**Hypothesis 2.6** implies that the transport term $f(\cdot; \mathcal{X})$ is nonnegative. It ensures the modeling assumptions stating that the phenotypic transfer is unidirectional: from healthy to cancer cells. The immediate consequence of this assumption is that the behavior of the system is mainly dictated by the location on the phenotype axis of the transport term support (i.e. the supports of $\pi$ and $\epsilon$).

Once the **Hypotheses 2.5 and 2.6** are assumed, the model dynamics can exhibit several types of behavior that can be associated to one of the following states:

* a stationary state where there is no cancer progression,

* abimodal state, where the cancer progression takes place but the system converges towards a bimodal distribution (i.e. an asymptotic state in which healthy cells and cancer cells are present),

* a pathological state where the healthy cells are totally replaced by cancer cells.

These three behaviors are determined by more refined assumptions on the functions involved in the system (4). The assumptions categorizing the dynamical behavior are summarized in Fig. 2.

**The stationary state.** The conditions ensuring this state are schematized in Fig. 2(a). Since there is no intersection between the support of the initial condition and the support of the transport term, there is no possible progression of the cancer. This behavior corresponds to the dynamics of the model when there are only healthy cells without proliferative cells (cf. Fig. 3(a)). Moreover, since there is no proliferative cells, the sensory axons density remains constant and the sympathetic axons density stays constant or follows a logistic growth (cf. Fig. 3(b)).

**The bimodal state.** The conditions ensuring the bimodal state are schematized in Fig. 2(b). The support of the initial density of cells on the phenotype axis is intersecting with the support of the functions implied in the progression towards the cancerous phenotype (the functions $\pi$ and $\eta$ in (1)). The function $\pi$ transports a proportion of the cell densities towards the cancer phenotype whereas the function $\eta$ starts accelerating the transport of the cell densities from the healthy state to the cancerous state only when cancer cells are present in the system. Nevertheless, in the bimodal state, the supports of $\pi$ and $\eta$ are not intersecting each other. In that case, the system dynamics tend towards a bimodal population distribution. The bimodal state corresponds to a state where there are cancer cells in the system after some time. However a proportion of non-cancerous cells is also present and persists through time (cf. Figs. 3(c) and 4(a)). Concerning the axons dynamics, in contrast to the stationary case, the sympathetic axons are first increasing thanks to the logistic growth and then decreasing thanks to the Allee effect (cf. Fig. 3(d)) and the sensory axons are following a logistic growth (since there are cancer cells in the system).

**The pathological state.** The asymptotic behavior corresponding to the pathological case is the opposite of the one from the stationary case. After a finite time, the cancer cells are present in the system and a progression towards cancerous phenotypes occurs for the remaining healthy cells (cf. Fig. 3(e)). Then, after some time, there are only cancer cells in the system. If a proportion of healthy cells remains after the start of the cancer progression (i.e. $\mathrm{supp}(Q_0) \cap \mathrm{supp}(\pi) \neq \emptyset$ and $\mathrm{supp}(\eta)/\mathrm{supp}(\pi) \neq \emptyset$), a second transport phenomenon towards the cancer phenotypes takes place because the supports of the function $\pi$ and the function $\eta$ are intersecting each other (cf. Figs. 2(c), 3(e) and 4(b)). This means that in this state, the cancer cells play a more prominent role in the cancer progression. Also, the cancer cell proliferation is controlled by the function $r$, its support on the phenotype axis and the various parameters involved in the saturation phenomena. In this case, since all cells are cancer cells after a finite time, the total amount
**Remark 3.** The Theorem 1 state the well-posedness on any finite time interval of the solution of (4). However, one can prove the global well-posedness with the same arguments assuming additional regularity assumptions (on the time variable $t \in \mathbb{R}_+$ instead of $t \in [0, T]$) on the functions $\theta$ and $r_{A_2}$ (e.g. no singularities). The main reason allowing the extension of the result about the well-posedness comes from the fact that assuming the Hypothesis 2.1 holds and $0 < A_1^0 < 1$ for $i = 1, 2$ then the trajectories $(\int_{\Omega} Q(\cdot, x)dx, A_1(t), A_2(t))$ remain bounded away from 0. However, finer stability estimates are required to conclude. Since this study focuses on the interactions between cells and axons during the tumorigenesis for the pancreatic cancer, the transient behaviour (instead of the asymptotic behaviour) is the main point of interest. The result of Theorem 1 is sufficient to pursue this goal.

## 2.3 Qualitative study of the system

The aim of this section is to study how the transport term (1) and the proliferation term (2) of the system (4) impact the final state of the model's dynamics. First, we show that, depending on the assumptions made about the functions $\pi$ (the basal speed of disease progression), $\eta$ (the speed of disease progression due to the presence of cancer cells), $r$ (the basal proliferation rate of cancer cells) and $Q_0$ (the initial distribution of cells with respect to their phenotypes), the different behaviours of the model's dynamics can be summed up in three categories. Secondly, we obtain explicit bounds on the first time of appearance of cells with cancerous phenotypes. This result is obtained under some assumptions on the transport term and relies on the characteristic curves of the equation.

### 2.3.1 Behaviour of the model

We recall that our phenotype domain $\Omega \subset \mathbb{R}$ can either be the whole real line or a segment centered in 0. Hence, a distinct separation for the cells with positive and negative phenotypes is assumed. The origin separates the cancer-induced proliferative cells from the healthy one. Moreover, at time $t = 0$ the cells are supposed to be at a non-cancerous and thus at a non-proliferating state. This assumption is formally translated in our model by the following hypothesis on the location on the phenotype axis of the initial distribution $Q_0$ (i.e. its support) and on the support of the proliferation rate $r$:

**Hypothesis 2.5.** Assume $Q_0$ is the initial condition of (4) which fulfills Hypothesis 2.1. Moreover, let $Q_0$ be compactly supported and

$$ \mathrm{supp}(Q_0) \cap \mathrm{supp}(r) = \emptyset. $$

Hypothesis 2.5 implies that the population of cells in a healthy pancreas remain constant until the apparition of cancer cells. Moreover, we also assume the following:

**Hypothesis 2.6.** Let

$$ \pi(x) \ge 0, \quad \eta(x) \ge 0, \quad \delta(x) \ge 0 \quad \text{and} \quad 1 - \beta(x)\rho(y) \ge 0 \quad \forall x \in \Omega, \forall y \in [0, 1]. $$

Hypothesis 2.6 implies that the transport term $f(\cdot; \mathcal{X})$ is nonnegative. It ensures the modeling assumptions stating that the phenotypic transfer is unidirectional: from healthy to cancer cells. The immediate consequence of this assumption is that the behaviour of the system is mainly dictated by the location on the phenotype axis of the transport term support (i.e. the supports of $\pi$ and $\epsilon$).

Once the Hypotheses 2.5 and 2.6 are assumed, the model dynamics can exhibit several types of behaviour that can be associated to one of the following states:

* a stationary state where there is no cancer progression,

* a bimodal state, where the cancer progression takes place but the system converges towards a bimodal distribution (i.e. an asymptotic state in which healthy cells and cancer cells are present),

* a pathological state where the healthy cells are totally replaced by cancer cells.

These three behaviours are determined by more refined assumptions on the functions involved in the system (4). The assumptions categorizing the dynamical behaviour are summarized in Figure 2.

**The stationary state.** The conditions ensuring this state are schematized in Figure 2a. Since there is no intersection between the support of the initial condition and the support of the transport term, there is no possible progression of the cancer. This behaviour corresponds to the dynamics of the model when there are only healthy cells without proliferative cells (cf. Figure 3a). Moreover, since there is no proliferative cells, the sensory axons density remains constant and the sympathetic axons density stays constant or follows a logistic growth (cf. Figure 3b).
<figure><img src="image_5.png" /><figcaption>(c) **The pathological state.** The support of the transport term (blue) is intersecting the support of the function $\eta$ (black) and the support of the initial condition $Q_0$ (red).</figcaption></figure>

<caption>Fig. 2. Scheme of the supports of the functions governing the dynamics of the cancer cells progression.</caption>

of cancer cells is higher as compared to the one in the bimodal state (cf. Figs. 3(c) and 3(e)). However, the axons dynamics in the bimodal state and the pathological state are similar since in both state cancer cells are present (cf. Fig. 3(f)).

### 2.3.2. Explicit bounds on the time to appearance of the first cancer cells

To obtain explicit bounds on the time to appearance of the first cancer cells, we state the following hypothesis which hold either on the bimodal case or the pathological case.

**Hypothesis 2.7.** Let

$$ \text{supp}(Q_0) \cap \text{supp}(\pi) \neq \emptyset, \quad \text{supp}(r) \cap \text{supp}(\pi) \neq \emptyset $$

and

$$ \pi: K \subseteq (\text{supp}(Q_0) \cap \text{supp}(\pi)) / \text{supp}(r) \to \mathbb{R}_+^*, \quad x \mapsto \pi(x) $$

be monotonous and increasing.

**Hypothesis 2.7** is stating that at least a proportion of cells from the initial condition progress towards the cancerous phenotype and that the cells can reach a proliferative state. Moreover, it also state that the function $\pi$ involved in the transport term is positive and increasing in a specific subset of its definition domain. The second statement ensures that the progression towards cancer cells happens in finite time and that progression is faster when the location on the phenotype axis is towards the pre-cancerous and cancerous states than when it is towards the healthy state. Hence, we consider the time $t^*$, the time to appearance of the first cancer cells, with the following definition.

**Definition 1.** We denote $t^*$ the first time when cells can proliferate, i.e.

$$ t^* := \inf \{ t \in \mathbb{R}_+^* | \min(\text{supp}(r)) \in \text{supp}(Q(t, \cdot)) \}. $$

Finally, we obtain an upper bound and a lower bound on the time $t^*$. The bounds depend on the assumptions made about the functions in the transport term ($\pi, \beta, \rho, \delta$ in (1)), those in the growth term ($r$ in (2)) and the initial condition $Q_0$ (cf. Figs. 2(b) or 2(c)).

**Proposition 1 (Appearance of the Proliferative Cells).** Assume Hypotheses 2.5, 2.6 and 2.7 hold. We denote (cf. Fig. 4)

$$ M_0 = \max \left( \text{supp}(Q_0) \cap \text{supp}(\pi) \right) \quad \text{and} \quad r_0 = \min(\text{supp}(r)). $$

Then there exist two positive constants $0 < m_1(\pi(M_0), \beta, \rho) < m_2(\pi, \delta)$ such that

$$ t^* \in \left[ \frac{r_0 - M_0}{m_2}; \frac{r_0 - M_0}{m_1} \right]. $$

**Proof (Proposition 1).** Let $t \in (0, t^*)$, then the following equation holds for $Q$:

$$ \partial_t Q(t, x) + \partial_x (f(t, x)Q(t, x)) = 0, $$

since $\text{supp}(Q_0) \cap \text{supp}(r) = \emptyset$. It implies that

$$ Q(t, x) = Q_0(X(0, t, x)) \times \exp \left( - \int_0^t \partial_x f \bigg|_{(t,x)=(s, X(s,t,x))} ds \right) $$

for $t \in (0, t^*)$ and $x \in \Omega$ where $X$ denote the characteristic curves, i.e.

$$ \left\{ \begin{array}{l} \frac{d}{ds} X(s, t, x) = f(s, X(s, t, x)), \\ X(t, t, x) = x. \end{array} \right. $$

Then, assuming $t \in (0, t^*)$, we have the following

$$ \text{supp}(Q(t, \cdot)) = \text{supp}(Q_0(X(0, t, \cdot))). $$

Hence, the bounds on $t^*$ can be found by studying the characteristics. Since **Hypothesis 2.6** hold, then there exist two functions $f_1: \Omega \to \mathbb{R}_+^*$ and $f_2: \Omega \to \mathbb{R}_+^*$ with

$$ f_1(y) = \pi(y)(1 - \beta(y)\rho(1)), $$

$$ f_2(y) = \pi(y)(1 + \delta(y)) $$

such that

$$ f_1(X(s, t, x)) \le \frac{d}{ds} X(s, t, x) \le f_2(X(s, t, x)). $$

In order to estimate the bounds on $t^*$, we focus on the characteristic starting at $M_0 = \max(\text{supp}(Q_0) \cap \text{supp}(\pi))$ at time 0 (cf. Fig. 4). Since $f(t, x) \ge 0$ for $t \ge 0$ and $x \in \Omega$ then

$$ X(t, 0, M_0) \ge M_0, \quad \text{for } t \ge 0. $$

On the first hand, there exists $\varepsilon > 0$ arbitrary small such that $\pi'(x) > 0$ for $x \in [M_0 - \varepsilon, M_0 + \varepsilon] \subset (\text{supp}(Q_0) \cap \text{supp}(\pi)) / \text{supp}(r)$ and $\max_{x \in \Omega} \beta(x)\rho(1) < 1$, we have that

$$ 0 < \pi(M_0) \left( 1 - \max_{x \in \Omega} \beta(x)\rho(1) \right) \le f_1(X(t, 0, M_0)), \quad \text{for } t \in (0, t^*). \quad (16) $$

On the other hand, since $t \in (0, t^*)$ then $N_c(t) = 0$, we have that

$$ f_2(X(t, 0, M_0)) \le \max_{x \in \Omega} \pi(x) \left( 1 + \max_{x \in \Omega} \delta(x) \right), \quad \text{for } t \in (0, t^*). \quad (17) $$

We denote $m_1 = \pi(M_0)(1 - \max_{x \in \Omega} \beta(x)\rho(1))$ and $m_2 = \max_{x \in \Omega} \pi(x)(1 + \max_{x \in \Omega} \delta(x))$, using the bounds in (16) and (17), we obtain

$$ M_0 + m_1 t \le X(t, 0, M_0) \le M_0 + m_2 t. $$

Finally, we obtain a bound for $t^*$ using the estimates on the characteristic starting at $M_0$:

$$ \frac{r_0 - M_0}{m_2} \le t^* \le \frac{r_0 - M_0}{m_1}. \quad \blacksquare $$
<figure><img src="image_4.png" /><figcaption>(c) The pathological state. The support of the transport term (blue) is intersecting the support of the function $\eta$ (black) and the support of the initial condition $Q_0$ (red).</figcaption></figure>

<caption>Figure 2: Scheme of the supports of the functions governing the dynamics of the cancer cells progression</caption>

**The bimodal state.** The conditions ensuring the bimodal state are schematized in Figure 2b. The support of the initial density of cells on the phenotype axis is intersecting with the support of the functions implied in the progression towards the cancerous phenotype (the functions $\pi$ and $\eta$ in (1)). The function $\pi$ transports a proportion of the cell densities toward the cancer phenotype whereas the function $\eta$ starts accelerating the transport of the cell densities from the healthy state to the cancerous state only when cancer cells are present in the system. Nevertheless, in the bimodal state, the supports of $\pi$ and $\eta$ are not intersecting each other. In that case, the system dynamics tend towards a bimodal population distribution. The bimodal state corresponds to a state where there are cancer cells in the system after some time. However a proportion of non-cancerous cells is also present and persists through time (cf. Figures 3c and 4a). Concerning the axons dynamics, in contrast to the stationary case, the sympathetic axons are first increasing thanks to the logistic growth and then decreasing thanks to the Allee effect (cf. Figure 3d) and the sensory axons are following a logistic growth (since there are cancer cells in the system).

**The pathological state.** The asymptotic behaviour corresponding to the pathological case is the opposite of the one from the stationary case. After a finite time, the cancer cells are present in the system and a progression towards cancerous phenotypes occurs for the remaining healthy cells (cf. Figure 3e). Then, after some time, there are only cancer cells in the system. If a proportion of healthy cells remains after the start of the cancer progression (i.e. $supp(Q_0) \cap supp(\pi) \neq \emptyset$ and $supp(\eta)/supp(\pi) \neq \emptyset$), a second transport phenomenon toward the cancer phenotypes takes place because the supports of the function $\pi$ and the function $\eta$ are intersecting each other (cf. Figures 2c, 3e and 4b). This means that in this state, the cancer cells play a more prominent role in the cancer progression. Also, the cancer cell proliferation is controlled by the function $r$, its support on the phenotype axis and the various parameters involved in the saturation phenomena. In this case, since all cells are cancer cells after a finite time, the total amount of cancer cells is higher as compared to the one in the bimodal state (cf. Figures 3c and 3e). However, the axons dynamics in the bimodal state and the pathological state are similar since in both state cancer cells are present (cf. Figure 3f).

### 2.3.2 Explicit bounds on the time to appearance of the first cancer cells

To obtain explicit bounds on the time to appearance of the first cancer cells, we state the following hypothesis which hold either on the bimodal case or the pathological case.

**Hypothesis 2.7.** Let

$$supp(Q_0) \cap supp(\pi) \neq \emptyset, \quad supp(r) \cap supp(\pi) \neq \emptyset$$

and

$$\pi: K \subsetneq (supp(Q_0) \cap supp(\pi))/supp(r) \to \mathbb{R}_+^*, \quad x \mapsto \pi(x)$$

be monotonous and increasing.

Hypothesis 2.7 is stating that at least a proportion of cells from the initial condition progress toward the cancerous phenotype and that the cells can reach a proliferative state. Moreover, it also states that the function $\pi$ involved in the transport term is positive and increasing in a specific subset of its definition domain. The second statement ensures that the progression toward cancer cells happens in finite time and that progression is faster when the location on the phenotype axis is towards the pre-cancerous and cancerous states than when it is towards the healthy state. Hence, we consider the time $t^*$, the time to appearance of the first cancer cells, with the following definition.

**Definition 1.** We denote $t^*$ the first time when cells can proliferate, i.e.

$$t^* := \inf \{t \in \mathbb{R}_+^* | \min(supp(r)) \in supp(Q(t, \cdot))\}.$$
<figure><img src="image_6.png" /><figcaption>(a) **The stationary state.** Evolution of the cell densities</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) **The bimodal state.** Evolution of the cell densities</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e) **The pathological state.** Evolution of the cell densities</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) **The stationary state.** Evolution of the axons</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) **The bimodal state.** Evolution of the axons</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(f) **The pathological state.** Evolution of the axons</figcaption></figure>

**Fig. 3.** Numerical simulations of the system (4) (cf. Appendix B for the details of the numerics). The sub Figs. 3(a), 3(c) and 3(f) show in the form of a heatmap the evolution of cell populations over time as a function of their phenotypes. The vertical axis is the time and the horizontal axis is the phenotype. The color map is an indicator of the cells densities. The sub Figs. 3(b), 3(d) and 3(f) describe the evolution of the normalized densities of sympathetic axons denoted $A_1$ and sensory axons denoted $A_2$.

The dynamics of the mathematical model illustrate the different states of pancreatic cancer tumorigenesis. Moreover, depending on the information available on the pro or anti-tumoral interactions between the axons and the cells, our model can be adapted to estimate the bounds on the time of appearance of the first cancer cells (Proposition 1).

### 3. Parametrization and modeling of the denervation

In order to confront our model to the data and be able to extract biological information from it, we need to parametrize the dynamical system and introduce degrees of freedom. The parametrization given in
<figure><img src="image_5.png" /><figcaption>(a) **The stationary state.** Evolution of the cell densities</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) **The stationary state.** Evolution of the axons</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) **The bimodal state.** Evolution of the cell densities</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) **The bimodal state.** Evolution of the axons</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(e) **The pathological state.** Evolution of the cell densities</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(f) **The pathological state.** Evolution of the axons</figcaption></figure>

Figure 3: Numerical simulations of the system (4) (cf. Appendix B for the details of the numerics). The subfigures 3a, 3c and 3f show in the form of a heatmap the evolution of cell populations over time as a function of their phenotypes. The vertical axis is the time and the horizontal axis is the phenotype. The color map is an indicator of the cells densities. The subfigures 3b, 3d and 3f describe the evolution of the normalized densities of sympathetic axons denoted $A_1$ and sensory axons denoted $A_2$.
<figure><img src="image_12.png" /><figcaption>**Fig. 4.** Scheme of the characteristic curves and the supports of the functions $\eta$, $Q_0$ and $\pi$ for the system (3) under the assumptions of Hypothesis 2.7. The red dashed lines are the characteristic curves. The red region is the region where the density of cells is non-zero. The thick red line is a dirac mass which contains part of the densities of the system (the characteristic curves are equal to 0). The time $t^*$ is the time of appearance of the first cancer cells (cf. Definition 1). The slope of the characteristic curves is changing at time $t^*$ because of the effect of the sensory axons and the effect of the cancer cells (through the function $\eta$) on the cancer progression.</figcaption></figure>

<table><caption>**Table 1**
Details on the parameters involved in the model.</caption><thead><tr><th></th><th>Name</th><th>Range/Value</th><th>Units</th><th>Description</th></tr></thead><tbody><tr><td rowspan="4">Discretization</td><td>T</td><td>70</td><td>day</td><td>life time of a mouse</td></tr><tr><td>L</td><td>50</td><td>-</td><td>phenotype domain</td></tr><tr><td>dt</td><td>min (<sup>0.9</sup>&frasl;<sub>&pi;<sub>0</sub>(1+&delta;)</sub>, 0.005)</td><td>day</td><td>time step</td></tr><tr><td>dx</td><td>0.05</td><td>-</td><td>space step</td></tr><tr><td rowspan="3">Initial condition</td><td><span style="text-decoration: overline;">Q</span><sub>0</sub></td><td>100</td><td>cells/mm<sup>3</sup></td><td>amount of healthy cells at time 0 days</td></tr><tr><td>A<sub>1</sub>(0)</td><td>0.15445</td><td>-</td><td>amount of sympathetic axons at time 0 days</td></tr><tr><td>A<sub>2</sub>(0)</td><td>0.004</td><td>-</td><td>amount of sensory axons at time 0 days</td></tr><tr><td rowspan="7">Speed of progression</td><td>&pi;<sub>0</sub></td><td>[0.1, 10]</td><td>day<sup>-1</sup></td><td>Basal amplitude of the speed phenotype at which progression starts</td></tr><tr><td>x<sub>1,x</sub></td><td>[30, 49]</td><td>-</td><td>steepness of the switch to progression</td></tr><tr><td>&epsilon;<sub>1,x</sub></td><td>[10<sup>-4</sup>, 10]</td><td>-</td><td>phenotype at which progression stops</td></tr><tr><td>x<sub>2,x</sub></td><td>10</td><td>-</td><td>steepness of the switch to no progression</td></tr><tr><td>&epsilon;<sub>2,x</sub></td><td>10</td><td>-</td><td>modulation by the sympathetic axons</td></tr><tr><td>&beta;</td><td>[10<sup>-4</sup>, 1]</td><td>-</td><td>modulation by the sensory axons</td></tr><tr><td>&delta;</td><td>[10<sup>-4</sup>, 1]</td><td>-</td><td>modulation by the sensory axons</td></tr><tr><td rowspan="5">Tumor growth</td><td>&gamma;<sub>r</sub></td><td>[10<sup>-2</sup>, 10]</td><td>day<sup>-1</sup></td><td>maximum proliferation rate</td></tr><tr><td>s<sub>r</sub></td><td>[0.2, 10]</td><td>-</td><td>rate of proliferation</td></tr><tr><td>&tau;<sub>c</sub></td><td>[120, 300]</td><td>cells/mm<sup>3</sup></td><td>carrying capacity of total concentration of cells</td></tr><tr><td>&mu;<sub>1</sub></td><td>[-1, 1]</td><td>-</td><td>modulation by the sympathetic axons</td></tr><tr><td>&mu;<sub>2</sub></td><td>[10<sup>-4</sup>, 1]</td><td>-</td><td>modulation by the sensory axons</td></tr><tr><td rowspan="5">Axon growth</td><td>A<sup>eq</sup><sub>1</sub></td><td>0.1544</td><td>-</td><td>Fixed from the data</td></tr><tr><td>r<sub>A<sub>1</sub></sub></td><td>[10<sup>-4</sup>, 5]</td><td>day<sup>-1</sup></td><td>Growth rate of the sympathetic axons</td></tr><tr><td>s<sub>&theta;</sub></td><td>[13, 20]</td><td>-</td><td>steepness of decrease</td></tr><tr><td>r&#773;<sub>A<sub>2</sub></sub></td><td>[10<sup>-4</sup>, 5]</td><td>day<sup>-1</sup></td><td>Growth rate of the sensory axons</td></tr><tr><td>s<sub>A<sub>2</sub></sub></td><td>[10<sup>-4</sup>, 10]</td><td>-</td><td>steepness of increase</td></tr></tbody></table>

Table 1 is described in Section 3.1. To understand the regulation coming from the nervous system, we also need to introduce a mathematical model of the denervation treatment. This is done in Section 3.2.

### 3.1. Details of the model and its parameters

**Choice of the domain.** We propose to use functions $\pi$ and $Q_0$ compactly supported so that we can reduce the study of the transport equation to a finite domain $\Omega = (-L, L)$ (cf. Section 2.3). We are then able to approximate the solution of the model using an upwind finite volume

scheme for the transport equation and explicit Euler scheme for the ODE (cf. Appendix B). In the following, we fix $L = 50$. We use a constant time step $dt$ and space step $h$ ensuring CFL conditions for the upwind scheme.

**Initial conditions.** We set $Q_0(x)$ to be the concentration of cells at time $t = 0$ days. Since initially, only healthy acinar cells are observed, and since negative values of $x$ close to 50 correspond to healthy cells, we then suppose $Q_0(x)$ follows a gaussian distribution centered at $-x_{init}$ defined by the following:
Finally, we obtain an upper bound and a lower bound on the time $t^*$. The bounds depend on the assumptions made about the functions in the transport term $(\pi, \beta, \rho, \delta$ in (1)), those in the growth term $(r$ in (2)) and the initial condition $Q_0$ (cf. Figures 2b or 2c).

<figure><img src="image_11.png" /><figcaption>Figure 4: Scheme of the characteristic curves and the supports of the functions $\eta, Q_0$ and $\pi$ for the system (3) under the assumptions of Hypothesis 2.7. The red dashed lines are the characteristic curves. The red region is the region where the density of cells is non-zero. The thick red line is a dirac mass which contains part of the densities of the system (the characteristic curves are equal to 0). The time $t^*$ is the time of appearance of the first cancer cells (cf. Definition 1). The slope of the characteristic curves is changing at time $t^*$ because of the effect of the sensory axons and the effect of the cancer cells (through the function $\eta$) on the cancer progression.</figcaption></figure>

**Proposition 1** (Appearance of the proliferative cells). Assume Hypotheses 2.5, 2.6 and 2.7 hold. We denote (cf. Figure 4)

$$M_0 = \max(\text{supp}(Q_0) \cap \text{supp}(\pi)) \quad \text{and} \quad r_0 = \min(\text{supp}(r)).$$

Then there exist two positive constants $0 < m_1(\pi(M_0), \beta, \rho) < m_2(\pi, \delta)$ such that

$$t^* \in \left[ \frac{r_0 - M_0}{m_2}, \frac{r_0 - M_0}{m_1} \right].$$

**Proof.**[Proposition 1] Let $t \in (0, t^*)$, then the following equation holds for $Q$:

$$\partial_t Q(t, x) + \partial_x (f(t, x)Q(t, x)) = 0,$$

since $\text{supp}(Q_0) \cap \text{supp}(r) = \emptyset$. It implies that

$$Q(t, x) = Q_0(X(0, t, x)) \times \exp \left( - \int_0^t \partial_x f \bigg|_{(t,x)=(s, X(s, t, x))} ds \right)$$

for $t \in (0, t^*)$ and $x \in \Omega$ where $X$ denote the characteristic curves, i.e.

$$\begin{cases} \frac{d}{ds} X(s, t, x) = f(s, X(s, t, x)), \\ X(t, t, x) = x. \end{cases}$$

Then, assuming $t \in (0, t^*)$, we have the following

$$\text{supp}(Q(t, \cdot)) = \text{supp}(Q_0(X(0, t, \cdot))).$$

Hence, the bounds on $t^*$ can be found by studying the characteristics. Since Hypothesis 2.6 hold, then there exist two functions $f_1 : \Omega \to \mathbb{R}_+^*$ and $f_2 : \Omega \to \mathbb{R}_+^*$ with

$$f_1(y) = \pi(y)(1 - \beta(y)\rho(1)),$$

$$f_2(y) = \pi(y)(1 + \delta(y))$$
<figure><img src="image_13.png" /><figcaption>Fig. 5. Illustration of the basal amplitude for the speed of disease progression $\pi$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 6. Illustration of the initial distribution of cells on the phenotype axis and the basal amplitude of the transport speed of the disease progression.</figcaption></figure>

$$Q_0(x) = \bar{Q}_0 \frac{\exp\left\{\frac{-(x+x_{init})^2}{4}\right\}}{\sqrt{4\pi}}$$

with $\bar{Q}_0 = 100$, and $x_{init} = 40$.

We note that $N(0)$ is close to $\bar{Q}_0$. Also, the choices of $\bar{Q}_0$, $L$ and $x_{init}$ are arbitrary and impact the values of many of the parameters.

As for the axons, at time $t = 0$ days, a negligible amount of sensory axons, and a small amount of autonomic axons are present. Moreover, the axon densities in the mathematical model (3) are renormalized and unitless. Based on the biological data provided, we consider

$$A_1(0) = 0.15445, \quad A_2(0) = 0.004.$$

**Cancer progression.** The experimental data that we have is not quantified on the impact that cancer cells have on cancer progression of the healthy cells in the late stages of the disease. Therefore, we impose in what follows $\eta = 0$ such that

$$f(t, x; \mathcal{X}) = \pi(x)(1 - \beta\rho(A_1) + \delta A_2).$$

This implies that transport to cancer phenotypes is no longer affected by the presence of cancer cells at a very advanced stage of the disease. The advantage of this modeling approach is that it reduces the number of parameters without significantly impacting the system's transient dynamics.

* The function $\pi$ represents a basal amplitude for the speed of disease progression:

$$\pi(x) = \pi_0 (\tanh(\epsilon_{1,\pi}(x + L - x_{1,\pi})) + \tanh(\epsilon_{2,\pi}(-x + L - x_{2,\pi}))) / 2.$$

Since initially the transdifferentiation of healthy cells is considered to be slow, we require the transport to be almost negligible around the lower bound of the phenotype domain. At some phenotype $-L + x_{1,\pi} \in (-L, 0)$, we expect to see an increase in the transport speed. A constant amplitude is then observed during the PanIN/PDAC phase. We impose the assumption that after $x = L - x_{2,\pi}$ we have decreasing transport speed in amplitude which vanishes in a neighborhood of the upper bound of the domain. This implies that no more transport takes place since cells have reached their PDAC phenotype for $x \ge L - x_{2,\pi}$ (cf. Fig. 5 and Fig. 6).

- Initially the amount of sensory axons is almost negligible and increases only after the presence of PDAC cancer cells. We then expect the speed before the presence of cancer cells to be

$$f(t, x; \mathcal{X}) \cong \pi_0(1 - \beta\rho(A_1)).$$

Since the cancer cells are seen around 35 days, the position of the initial distribution on the phenotype axis provides information on the approximate value of the parameters linked to cancer progression and the following can be considered

$$\pi_0(1 - \beta\rho(A_1)) \cong \frac{x_{init}}{35}.$$

But as $\rho(A_1) < 1$, we then have

$$\pi_0(1 - \beta) \le \pi_0(1 - \beta\rho(A_1)) \le \pi_0$$

which gives,

$$\pi_0 \in \left[ \frac{x_{init}}{35}, \frac{x_{init}}{35(1 - \beta)} \right].$$

- The parameter $x_{1,\pi}$ is closely related to the location on the phenotype axis where the cancer progression starts. We note that based on Section 2.3, we would like to secure that $\text{supp}(Q_0) \cap \text{supp}(\pi) \ne \emptyset$ so that tumorigenesis takes place. Since $Q_0$ is constructed such that $Q_0(x) \approx 0$ for $x \ge -20$, we would then like to ensure that $\min\{\text{supp}(\pi)\} \le -30$. Thus, we choose $L - x_{1,\pi} \in [30, 49]$ so that $x_{1,\pi} \in [1, 20]$.

- The parameter $\epsilon_{1,\pi}$ determines the rapidity of increase in the speed of transport. No biological information has been provided on this value but keeping in mind that it is essential for $\text{supp}(Q_0) \cap \text{supp}(\pi) \ne \emptyset$, we then take $\epsilon_{1,\pi} \in [10^{-4}, 10]$.

- The parameter $x_{2,\pi}$ is linked to the location on the phenotype axis where the transport decreases and stops. The idea of introducing $x_{2,\pi}$ is to avoid losing mass at the boundary value $L$. Thus, we choose $x_{2,\pi}$ so that $L - x_{2,\pi} < 50$. Here, we fix this value to $x_{2,\pi} = 10$.

- The parameter $\epsilon_{2,\pi}$ determines the speed at which the transport speed is reduced. Since we do not expect any transport to happen after all susceptible cells have attained their cancerous state, we then choose $\epsilon_{2,\pi} = 10$ so that $\pi(x)$ vanishes fast after $L - x_{2,\pi}$.

* The parameters $\beta$ and $\delta$ represent the maximum rate of the regulation of the progression by the sympathetic and sensory axons respectively. Since no biological information is provided about their value and since they represent a modulation, we choose $(\beta, \delta) \in [10^{-4}, 1] \times [10^{-4}, 1]$.

* The function $\rho$ is a function that removes the effect of regulation for small values of $A_1$ and that is linear for larger values of $A_1$. This regulation through the function $\rho$ is motivated by the fact that the autonomous axons are non-negligibly present in healthy pancreas. We denote $A_1^{eq}$ the renormalized value of $A_1$ at the healthy state. To ensure that the quantity $A_1$ does not artificially affect the cancer progression when $A_1 < A_1^{eq}$, we construct $\rho$ as the following:

<figure><img src="image_15.png" /></figure>

**Cancer growth.** The proliferation is generated by a logistic growth. More precisely,

$$g(t, x; \mathcal{X}) := r(x)Q(t, x) \left( 1 - \frac{N(t)}{\tau_C} - \mu_1 A_1(t) + \mu_2 A_2(t) \right).$$
such that

$$f_1(X(s, t, x)) \leq \frac{d}{ds} X(s, t, x) \leq f_2(X(s, t, x)).$$

In order to estimate the bounds on $t^*$, we focus on the characteristic starting at $M_0 = \max(\text{supp}(Q_0) \cap \text{supp}(\pi))$ at time 0 (cf. Figure 4). Since $f(t, x) \geq 0$ for $t \geq 0$ and $x \in \Omega$ then

$$X(t, 0, M_0) \geq M_0, \quad \text{for } t \geq 0.$$

On the first hand, there exists $\varepsilon > 0$ arbitrary small such that $\pi'(x) > 0$ for $x \in [M_0 - \varepsilon, M_0 + \varepsilon] \subset (\text{supp}(Q_0) \cap \text{supp}(\pi)) / \text{supp}(r)$ and $\max_{x \in \Omega} \beta(x)\rho(1) < 1$, we have that

$$0 < \pi(M_0)(1 - \max_{x \in \Omega} \beta(x)\rho(1)) \leq f_1(X(t, 0, M_0)), \quad \text{for } t \in (0, t^*). \qquad (9)$$

On the other hand, since $t \in (0, t^*)$ then $N_c(t) = 0$, we have that

$$f_2(X(t, 0, M_0) \leq \max_{x \in \Omega} \pi(x)(1 + \max_{x \in \Omega} \delta(x)), \quad \text{for } t \in (0, t^*). \qquad (10)$$

We denote $m_1 = \pi(M_0)(1 - \max_{x \in \Omega} \beta(x)\rho(1))$ and $m_2 = \max_{x \in \Omega} \pi(x)(1 + \max_{x \in \Omega} \delta(x))$, using the bounds in (9) and (10), we obtain

$$M_0 + m_1 t \leq X(t, 0, M_0) \leq M_0 + m_2 t.$$

Finally, we obtain a bound for $t^*$ using the estimates on the characteristic starting at $M_0$:

$$\frac{r_0 - M_0}{m_2} \leq t^* \leq \frac{r_0 - M_0}{m_1}.$$

<figure><img src="image_12.png" /></figure>

The dynamics of the mathematical model illustrate the different states of pancreatic cancer tumorigenesis. Moreover, depending on the information available on the pro or anti-tumoral interactions between the axons and the cells, our model can be adapted to estimate the bounds on the time of appearance of the first cancer cells (Proposition 1).

# 3 Parametrization and modelling of the denervation

In order to confront our model to the data and be able to extract biological information from it, we need to parametrize the dynamical system and introduce degrees of freedom. The parametrization given in Table 1 is described in Section 3.1. To understand the regulation coming from the nervous system, we also need to introduce a mathematical model of the denervation treatment. This is done in Section 3.2.

## 3.1 Details of the model and its parameters

**Choice of the domain.** We propose to use functions $\pi$ and $Q_0$ compactly supported so that we can reduce the study of the transport equation to a finite domain $\Omega = (-L, L)$ (cf. Section 2.3). We are then able to approximate the solution of the model using an upwind finite volume scheme for the transport equation and explicit Euler scheme for the ODE (cf. Appendix B). In the following, we fix $L = 50$. We use a constant time step $dt$ and space step $h$ ensuring CFL conditions for the upwind scheme.

**Initial Conditions.** We set $Q_0(x)$ to be the concentration of cells at time $t = 0$ days. Since initially, only healthy acinar cells are observed, and since negative values of $x$ close to 50 correspond to healthy cells, we then suppose $Q_0(x)$ follows a gaussian distribution centered at $-x_{init}$ defined by the following:

$$Q_0(x) = \bar{Q}_0 \frac{\exp\left\{\frac{-(x+x_{init})^2}{4}\right\}}{\sqrt{4\pi}}$$

with $\bar{Q}_0 = 100$, and $x_{init} = 40$.

We note that $N(0)$ is close to $\bar{Q}_0$. Also, the choices of $\bar{Q}_0$, $L$ and $x_{init}$ are arbitrary and impact the values of many of the parameters.

As for the axons, at time $t = 0$ days, a negligible amount of sensory axons, and a small amount of autonomic axons are present. Moreover, the axon densities in the mathematical model (3) are renormalized and unitless. Based on the biological data provided, we consider

$$A_1(0) = 0.15445, \quad A_2(0) = 0.004.$$
<figure><img src="image_16.png" /><figcaption>Fig. 7. Illustration of the behavior of the basal growth rate of the pre-tumor and tumor cells.</figcaption></figure>

* The function $r$ represents the basal growth rate of the pre-tumor and tumor cells (see Fig. 7). We expect for proliferation to start taking place in the presence of pre-cancerous cells, i.e. in a neighborhood of 0$^{-}$. Hence, we consider that

$$r(x) = \frac{1}{2} \gamma_r (1 + \tanh(s_r x)).$$

- The maximum proliferation rate that the cells can exhibit is $\gamma_r$. Here, we take $\gamma_r \in [10^{-2}, 10]$.

- The parameter $s_r$ determines the rate at which proliferation takes place. For small values of $s_r$, then the proliferation increases in rate in precancerous lesions, whereas for big values of $s_r$, the tumor growth increases in rate only in late stages of precancerous lesions. We choose $s_r \in [0.2, 10]$.

* The growth rate of the pre-tumor and tumor cells can be regulated by the presence of axons through the term

$$-\mu_1 A_1(t) + \mu_2 A_2(t).$$

In that, we consider that the sympathetic axons either slow down or accelerate tumor growth whereas the sensory axons only accelerate tumor growth. The parameters $\mu_1$ and $\mu_2$ modulate the impact that the axons have on the proliferation. We consider that $\mu_1 \in [-1, 1]$ and $\mu_2 \in [10^{-4}, 1]$.

* It is important to point out that even if cells proliferate, the total quantity of cells is limited and there is no explosion in finite time. In other words, we require $1 - \frac{N(t)}{\tau_C} - \mu_1 A_1(t) + \mu_2 A_2(t) < \infty$ with $\tau_C$ being the carrying capacity of the amount of cells in the absence of axons. Thus, using the fact that $A_1 \le 1$ for $i = 1, 2$ we have,

$$\begin{aligned} N(t) &\le \tau_C(1 - \mu_1 A_1(t) + \mu_2 A_2(t)) \le \tau_C(1 + \mu_2 + \mu_1^{-}), \\ &\le C(\tau_C), \end{aligned}$$

with $C(\tau_C) = \tau_C(1 + \mu_2 + \mu_1^{-})$ where $x^{-} = \frac{1}{2}(|x| - x)$. Since $N(0) = 100$, we then require $N(0) \le \tau_C$. We choose $\tau_C \in [120, 300]$.

**Sympathetic axons growth.** We introduced an Allee effect in the growth term of the sympathetic axons (similar to what is seen in Lolas et al., 2016). We recall that the dynamics of $A_1$ is given by the following

$$r_{A_1} A_1(t) \left( \frac{A_1(t)}{\theta \left( \frac{N(t)}{N(0)} \right)} - 1 \right) (1 - A_1(t)).$$

* The parameter $r_{A_1}$ is the growth rate of the autonomic axons. The bigger the value of $r_{A_1}$ is, the steeper the increase or decrease of $A_1$ is. We consider $r_{A_1} \in [10^{-4}, 5]$.

* In the light of what is done in Lolas et al. (2016), we also propose to modulate the threshold $\theta$ with the total amount of cells:

$$\theta \left( \frac{N(t)}{N(0)} \right) = \frac{A_1^{eq}}{2} + \frac{1}{2} \tanh \left( s_0 \left( \frac{N(t)}{N(0)} - 1 - 0.1 \right) \right) + \frac{1}{2} \quad (18)$$

where $A_1^{eq}$ is the density of the sympathetic axons at the healthy state and $s_0$ is the steepness of the decrease speed of $A_1$. The value of $A_1^{eq}$ is estimated from the data (cf. Table 1). In order to ensure a biological meaning, the following has to hold: for small time $t$, $\theta \left( \frac{N(t)}{N(0)} \right) \approx \theta(1)$,

and we want $\theta(1) < A_1(0)$ with $A_1(0) \approx A_1^{eq}$ to observe initially an increase in the autonomic axon density. Thus, we want

$$s_{\theta} \ge \frac{\tanh^{-1}(1 - A_1^{eq})}{10^{-1}}.$$

It leads us to consider $s_\theta \in [13, 20]$. For large values of the ratio $\frac{N(t)}{N(0)}$, we have $\theta \left( \frac{N(t)}{N(0)} \right) \approx 1 + \frac{A_1^{eq}}{2} > 1$ leading to a decrease in the autonomic axon density.

**Sensory axons growth.** The sensory axons follow a logistic growth regulated by the presence of PDAC

$$r_{A_2} \left( \frac{N_c(t)}{N(t)} \right) A_2(1 - A_2)$$

with $r_{A_2}(s) = \bar{r}_{A_2} \tanh(s_{A_2}s)$. The parameter $\bar{r}_{A_2}$ is the maximum growth rate of the sensory axons. We consider $\bar{r}_{A_2} \in [10^{-4}, 5]$. The parameter $s_{A_2}$ modulates the steepness of the increase of the growth of $A_2$. Here, we consider $s_{A_2} \in [10^{-4}, 10]$.

## 3.2. Modeling denervation treatment

Testing the effects of denervation in preclinical models often requires the use of numerous animals, raising ethical concerns. Therefore, the use of a mathematical model to predict the effects of denervation becomes highly valuable, allowing for the replacement and reduction of the number of animals used. In the literature, the role of axons in cancer progression has been studied through chemical or surgical denervation. In our mathematical model, denervation is implemented by setting parameters related to the targeted axon to 0 after the time of the intervention.

Specifically, denervating sympathetic axons is equivalent to suppressing the influence of $A_1$ on cell transport and proliferation. This is achieved by setting the parameters as follows after the denervation time:

$$\beta = 0, \quad \mu_1 = 0.$$

Conversely, denervating sensory axons is equivalent to suppressing the influence of $A_2$ on cell transport and proliferation, after the denervation time, the parameters are set as:

$$\delta = 0, \quad \mu_2 = 0.$$

For denervating both axons, the parameters are set as follows after the denervation time:

$$\beta = 0, \quad \mu_1 = 0, \quad \delta = 0, \quad \mu_2 = 0.$$

We have to adapt the notations given above to take into account the above mathematical denervation when useful. The dynamical system defined in Section 2, in particular in Eq. (3) is parameterized by $\theta$ as explained in Section 3.1 and models the control treatment (AA). The total amount of cells and the total amount of cancer cells arising from this dynamical system is denoted by $N(t|\theta, AA)$ and $N_c(t|\theta, AA)$, respectively. The dynamical system that models the denervation treatment is obtained from the same dynamical system, yet the value of $\theta$ is changed after time of the intervention, where some components are set to 0 as explained above. The total amount of cells and the total amount of cancer cells arising from this second dynamical system is denoted by $N(t|\theta, OHDA)$ and $N_c(t|\theta, OHDA)$, respectively.

## 4. Calibration

In this section, we focus on the calibration of the model to the biological data at our disposal. It is rather complex in our case since

* a small number of biological data is provided, with variability in samples as described in Section 4.1,
<table><caption>Table 1: Details on the parameters involved in the model.</caption><thead><tr><th></th><th>Name</th><th>Range/Value</th><th>Units</th><th>Description</th></tr></thead><tbody><tr><td rowspan="4">Discretization</td><td>T</td><td>70</td><td>day</td><td>life time of a mouse</td></tr><tr><td>L</td><td>50</td><td>-</td><td>phenotype domain</td></tr><tr><td>dt</td><td>min(<sup>0.9dx</sup>&frasl;<sub>&pi;<sub>0</sub>(1+&delta;)</sub>, 0.005)</td><td>day</td><td>time step</td></tr><tr><td>dx</td><td>0.05</td><td>-</td><td>space step</td></tr><tr><td rowspan="3">Initial condition</td><td>Q<sub>0</sub></td><td>100</td><td>cells/mm<sup>3</sup></td><td>amount of healthy cells at time 0 days</td></tr><tr><td>A<sub>1</sub>(0)</td><td>0.15445</td><td>-</td><td>amount of sympathetic axons at time 0 days</td></tr><tr><td>A<sub>2</sub>(0)</td><td>0.004</td><td>-</td><td>amount of sensory axons at time 0 days</td></tr><tr><td rowspan="7">Speed of progression</td><td>&pi;<sub>0</sub></td><td>[0.1, 10]</td><td>day<sup>-1</sup></td><td>Basal amplitude of the speed phenotype at which progression starts</td></tr><tr><td>x<sub>1,&pi;</sub></td><td>[30, 49]</td><td>-</td><td>steepness of the switch to progression</td></tr><tr><td>&epsilon;<sub>1,&pi;</sub></td><td>[10<sup>-4</sup>, 10]</td><td>-</td><td>phenotype at which progression stops</td></tr><tr><td>x<sub>2,&pi;</sub></td><td>10</td><td>-</td><td>steepness of the switch to no progression</td></tr><tr><td>&epsilon;<sub>2,&pi;</sub></td><td>10</td><td>-</td><td>modulation by the sympathetic axons</td></tr><tr><td>&beta;</td><td>[10<sup>-4</sup>, 1]</td><td>-</td><td>modulation by the sensory axons</td></tr><tr><td>&delta;</td><td>[10<sup>-4</sup>, 1]</td><td>-</td><td></td></tr><tr><td rowspan="5">Tumor growth</td><td>&gamma;<sub>r</sub></td><td>[10<sup>-2</sup>, 10]</td><td>day<sup>-1</sup></td><td>maximum proliferation rate</td></tr><tr><td>s<sub>r</sub></td><td>[0.2, 10]</td><td>-</td><td>rate of proliferation</td></tr><tr><td>&tau;<sub>c</sub></td><td>[120, 300]</td><td>cells/mm<sup>3</sup></td><td>carrying capacity of total concentration of cells</td></tr><tr><td>&mu;<sub>1</sub></td><td>[-1, 1]</td><td>-</td><td>modulation by the sympathetic axons</td></tr><tr><td>&mu;<sub>2</sub></td><td>[10<sup>-4</sup>, 1]</td><td>-</td><td>modulation by the sensory axons</td></tr><tr><td rowspan="5">Axon growth</td><td>A<sup>eq</sup><sub>1</sub></td><td>0.1544</td><td>-</td><td>Fixed from the data</td></tr><tr><td>r<sub>A<sub>1</sub></sub></td><td>[10<sup>-4</sup>, 5]</td><td>day<sup>-1</sup></td><td>Growth rate of the sympathetic axons</td></tr><tr><td>s<sub>&theta;</sub></td><td>[13, 20]</td><td>-</td><td>steepness of decrease</td></tr><tr><td>r&#773;<sub>A<sub>2</sub></sub></td><td>[10<sup>-4</sup>, 5]</td><td>day<sup>-1</sup></td><td>Growth rate of the sensory axons</td></tr><tr><td>s<sub>A<sub>2</sub></sub></td><td>[10<sup>-4</sup>, 10]</td><td>-</td><td>steepness of increase</td></tr></tbody></table>

**Cancer progression.** The experimental data that we have is not quantified on the impact that cancer cells have on cancer progression of the healthy cells in the late stages of the disease. Therefore, we impose in what follows $\eta = 0$ such that

$$f(t, x; \mathcal{X}) = \pi(x)(1 - \beta\rho(A_1) + \delta A_2).$$

This implies that transport to cancer phenotypes is no longer affected by the presence of cancer cells at a very advanced stage of the disease. The advantage of this modeling approach is that it reduces the number of parameters without significantly impacting the system's transient dynamics.

<figure><img src="image_13.png" /><figcaption>Figure 5: Illustration of the basal amplitude for the speed of disease progression $\pi$.</figcaption></figure>

* The function $\pi$ represents a basal amplitude for the speed of disease progression:

$$\pi(x) = \pi_0 (\tanh(\epsilon_{1,\pi} (x + L - x_{1,\pi})) + \tanh(\epsilon_{2,\pi} (-x + L - x_{2,\pi}))) / 2.$$

Since initially the transdifferentiation of healthy cells is considered to be slow, we require the transport to be almost negligible around the lower bound of the phenotype domain. At some phenotype $-L + x_{1,\pi} \in (-L, 0)$, we expect to see an increase in the transport speed. A constant amplitude is then observed during the PanIN/PDAC phase. We impose the assumption that after $x = L - x_{2,\pi}$ we have decreasing transport speed in amplitude which vanishes in a neighborhood of the upper bound of the domain. This implies that no more transport takes place since cells have reached their PDAC phenotype for $x \ge L - x_{2,\pi}$ (cf. Figure 5).

- Initially the amount of sensory axons is almost negligible and increases only after the presence of PDAC cancer cells. We then expect the speed before the presence of cancer cells to be

$$f(t, x; \mathcal{X}) \cong \pi_0(1 - \beta\rho(A_1)).$$
* the model is far from being linear in the parameters and the predictions returned by the mathematical model cannot be written explicitly as a function of the parameters,

* we do not have at our disposal a likelihood function to explain the difference between the predictions of model and the data.

Moreover, the set of parameters is of rather high dimension ($d = 14$) and we need to avoid selecting the best set of parameters that would overfit the few observed data and would not reflect the biological phenomenon we are trying to capture.

To this aim, we propose in Section 4.2 a 2-dimension criterion to measure the quality of the calibration of the model to the data and to our biological knowledge explained in Section 4.1: (1) a mean squared error criterion based on the cell distribution and axon evolution, that measures the difference between model outcomes and data, and (2) a criterion based on the biological knowledge we have on the chronology of cell evolution. As described in Section 4.3, we managed to extract biologically relevant sets of parameters that are consistent with the cell data and the chronological knowledge by exploring the parameter space with a multi-stage quasi-Monte Carlo algorithm that starts from a flat, uniform prior over the parameter space which is the hypercube given by the range column of Table 1, and ends with the 0.4% best sets of parameters according to the criterion, that are presented in Section 4.4 and in Appendix B.2.

## 4.1. Biological data

The selection of experimental data for calibrating the mathematical model is crucial and should align with the characteristics of the relevant preclinical model. The calibration process encompasses various aspects of the animal model, incorporating chronological insights into cancer progression derived from the literature. Additionally, data on sympathetic nerve axons density, sensory nerve axons density, and cell density at different stages of cancer progression are considered during the calibration process.

**Mice model and data acquisition methods.** In this study, we use the *Kras*<sup>LSL-G12D/+</sup>; *Cdkn2a*(*Ink4a/Arf*)<sup>lox/lox</sup>; *Pdx1-Cre* (KIC) mouse model of pancreatic cancer. The overexpression of a mutated form of the oncogene *KRAS* in pancreatic cells induces transdifferentiation, where healthy acinar cells (the functional unit of the exocrine pancreas) transform into ductal cell-like cells, leading to the formation of acinar to ductal metaplasia (ADM). Subsequently, ADM can progress to form premalignant pancreatic intraepithelial neoplasia (PanIN). The deletion of tumor-suppressor genes *Ink4 A* and *Arf* further accelerates tumor development, resulting in PDAC itself. The KIC mice model successfully replicates the stepwise progression observed in human pancreatic cancer (cf. *Aguirre et al.,* 2003, and Fig. 8). Several mouse pancreases at the time points of 6.5 weeks and 8 weeks are utilized in the experiments to accommodate biological heterogeneity. The technique used to quantify nerve axons density and cell density within lobules of the pancreatic tissue is detailed in the source data file of *Guillot et al.* (2022). In short, data comes from the quantification of 3D cleared tissue images obtained by Light Sheet Fluorescence Microscopy. Tissue staining was achieved through immunostaining and tissue clearing using the iDISCO+ protocol. Tyrosine hydroxylase antibody staining (TH), an enzyme involved in the biosynthesis of norepinephrine (one of the main neurotransmitters of the sympathetic nervous system), and Calcitonin gene-related peptide (CGRP) antibody were used to visualize sympathetic and sensory innervation respectively. Regions of Interest (ROI), namely asymptomatic (ASYM), ADM, PanIN, and PDAC, were segmented based on the autofluorescence signal of the tissue, and their volumes were measured. Axonal networks were manually reconstructed to measure the total axon length in each ROI. The unitless axon density was calculated as follows: $$ \left( \frac{\text{axons length sum} \times 1000}{\text{µm}} \right) / \text{(volume of each Regions of Interest (ROI))}^{\frac{1}{3}} (\text{µm}) $$

. Densities obtained through experimentation were normalized by first dividing by the cubic root of the volume in which the density was measured and then by an affine transformation to obtain proportions. This normalization facilitates the comparison of quantities both within the experimental data and with the outputs provided by the mathematical model (cf. Figs. 9 and 10).

The method used to describe the proportion of tissue categorized by their phenotypes in histological sections through the pancreas of a 6.5-week-old mouse can be found in the source data file of *Guillot et al.* (2022) (cf. Fig. 11).

**Knowledge on the chronological process.** Based on the previous characterization of the KIC model (seen in *Aguirre et al.,* 2003 and in our personal observation), the time of observation of ADM lobules denoted $t_{ADM}$ is around 17 days, the first appearance of relatively advanced-staged PanIN lobules denoted $t_{PANIN}$ is observed around 24.5 days, and the first appearance of PDAC lobules denoted $t_{PDAC}^{early}$ is observed around 35 days. The median survival of a KIC mouse model is 9 weeks and so, at 6.5 weeks ($= t_{PDAC}$) the mice has already developed a tumor, but at 8 weeks ($= t_{PDAC}^{advanced}$), the cancer is really aggressive. All these chronological assumptions are summarized in Fig. 8. Additionally, we include a time interval of plus or minus three days around the above temporal values to account for variability in observations.

**Sympathetic and sensory nerve axons quantification.** The normalized densities of sympathetic and sensory nerve axons are shown in Fig. 9 and Fig. 10, respectively, quantified at 6.5 weeks of the KIC age. For sensory axons, density measurements were also performed at 8 weeks for the PDAC lesions. Also, each sample of density measurements from the data does not only characterize the axon densities but also the tissue stages from where it was observed (ROI). This additional information may not be fully exploited by the model because information about space is not taken into account in the construction of the model. One way to use all information of the data is to extrapolate a dynamic based on the ROI type of the data origin. Since we know that $t_{ASYM}$, $t_{ADM}$, $t_{PANIN}$, and $t_{PDAC}$ are the times of observation of each ROI type during cancer progression, we then consider that the quantification of the normalized densities of the axons in each ROI type corresponds to those times respectively. Those timings are stated on the x-axis of Figs. 9 and 10. Hence, the dynamics of the axons are illustrated by the green curves in Figs. 9 and 10. It allows to combine qualitative information (chronological knowledge) to quantitative information and ensures a more precise calibration of the mathematical model.

**PDAC cell concentration.** Fig. 11 illustrates the data extracted for the normalized cancer cell concentration at 6.5 weeks. We set the time of observation for the data and aim to identify model dynamics that closely match these data points at 6.5 weeks. Through observation, we note that the proportion of PDAC cells is higher in OHDA samples than in AA samples.

## 4.2. A 2-dimensional criterion

The first dimension of our criterion is a mean squared error criterion that measures the difference between model outcomes and data. We consider the following data:

* n biological samples of cell distributions $p_i \in (0; 1)$ under treatment $u_i \in \{AA, OHDA\}$ ($i = 1, \dots, n$) that should be compared to the proportion of proliferating cells $p(t|\theta, u)$ given by the model, on average over time $t \in (0, T)$

* m biological samples of axon density $a_i \in (0; 1)$ in AA ($i = n + 1, \dots, n + m$) that should be compared to the axon density evolution $A(t|\theta)$ given by the model, on average over the time range $t \in I_i$ corresponding to the kind of tissues of *ith* biological sample (ADM, PanIN, PDAC or mature PDAC). Moreover, $z_i \in \{A_1, A_2\}$ indicates whether the observed density $a_i$ is a density of sympathetic ($A_1$) or sensory ($A_2$) axons.
Since the cancer cells are seen around 35 days, the position of the initial distribution on the phenotype axis provides information on the approximate value of the parameters linked to cancer progression and the following can be considered

$$ \pi_0(1 - \beta\rho(A_1)) \approx \frac{x_{init}}{35}. $$

But as $\rho(A_1) < 1$, we then have

$$ \pi_0(1 - \beta) \le \pi_0(1 - \beta\rho(A_1)) \le \pi_0 $$

which gives,

$$ \pi_0 \in \left[ \frac{x_{init}}{35}, \frac{x_{init}}{35(1 - \beta)} \right]. $$

- The parameter $x_{1,\pi}$ is closely related to the location on the phenotype axis where the cancer progression starts. We note that based on Section 2.3, we would like to secure that $\text{supp}(Q_0) \cap \text{supp}(\pi) \neq \emptyset$ so that tumorigenesis takes place. Since $Q_0$ is constructed such that $Q_0(x) \approx 0$ for $x \ge -20$, we would then like to ensure that $\min\{\text{supp}(\pi)\} \le -30$. Thus, we choose $L - x_{1,\pi} \in [30, 49]$ so that $x_{1,\pi} \in [1, 20]$.

<figure><img src="image_14.png" /><figcaption>Figure 6: Illustration of the initial distribution of cells on the phenotype axis and the basal amplitude of the transport speed of the disease progression.</figcaption></figure>

- The parameter $\epsilon_{1,\pi}$ determines the rapidity of increase in the speed of transport. No biological information has been provided on this value but keeping in mind that it is essential for $\text{supp}(Q_0) \cap \text{supp}(\pi) \neq \emptyset$, we then take $\epsilon_{1,\pi} \in [10^{-4}, 10]$.

- The parameter $x_{2,\pi}$ is linked to the location on the phenotype axis where the transport decreases and stops. The idea of introducing $x_{2,\pi}$ is to avoid losing mass at the boundary value $L$. Thus, we choose $x_{2,\pi}$ so that $L - x_{2,\pi} < 50$. Here, we fix this value to $x_{2,\pi} = 10$.

- The parameter $\epsilon_{2,\pi}$ determines the speed at which the transport speed is reduced. Since we do not expect any transport to happen after all susceptible cells have attained their cancerous state, we then choose $\epsilon_{2,\pi} = 10$ so that $\pi(x)$ vanishes fast after $L - x_{2,\pi}$.

* The parameters $\beta$ and $\delta$ represent the maximum rate of the regulation of the progression by the sympathetic and sensory axons respectively. Since no biological information is provided about their value and since they represent a modulation, we choose $(\beta, \delta) \in [10^{-4}, 1] \times [10^{-4}, 1]$.

* The function $\rho$ is a function that removes the effect of regulation for small values of $A_1$ and that is linear for larger values of $A_1$. This regulation through the function $\rho$ is motivated by the fact that the autonomous axons are non-negligibly present in healthy pancreas. We denote $A_1^{eq}$ the renormalized value of $A_1$ at the healthy state. To ensure that the quantity $A_1$ does not artificially affect the cancer progression when $A_1 < A_1^{eq}$, we construct $\rho$ as the following:

$$ \rho(x) = \begin{cases} 0 & \text{for } 0 \le x \le A_1^{eq}, \\ y = x - A_1^{eq} & \text{for } A_1^{eq} \le x \le 1. \end{cases} $$

<figure><img src="image_15.png" /></figure>

**Cancer growth.** The proliferation is generated by a logistic growth. More precisely,

$$ g(t, x; \mathcal{X}) := r(x)Q(t, x) \left( 1 - \frac{N(t)}{\tau_C} - \mu_1 A_1(t) + \mu_2 A_2(t) \right). $$

* The function $r$ represents the basal growth rate of the pre-tumor and tumor cells. We expect for proliferation to start taking place in the presence of pre-cancerous cells, i.e. in a neighborhood of $0^{-}$. Hence, we consider that

$$ r(x) = \frac{1}{2}\gamma_r(1 + \tanh(s_r x)). $$

- The maximum proliferation rate that the cells can exhibit is $\gamma_r$. Here, we take $\gamma_r \in [10^{-2}, 10]$.
<figure><img src="image_17.png" /><figcaption>**Fig. 8.** Diagram outlining the chronological assumptions regarding cancer progression in pancreas of KIC mice, drawing from Aguirre et al. (2003) and our firsthand observations. The indicated times represent the estimated average time at which each new stage is expected to be observed for the first time. The bold red segments indicate a time interval of plus or minus three days around the mentioned temporal values, encompassing variations in chronological estimates.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 9.** Whisker Plot of Normalized Sympathetic Nerve Axon Density: Data from different KIC mice depict the sympathetic axonal density within various ROI in the pancreas at 6.5 weeks of age, categorized as ASYMP, ADM, PanIN, or PDAC. The quantification involved 4 ASYMP, 5 ADM, 6 PanIN, and 7 PDAC ROIs (Guillot et al., 2022). The dotted blue points correspond to the data. The green curve represents the mean density in each ROI, with identified stages on the x-axis extrapolated from the chronological knowledge of each ROI type.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 10.** Whisker Plot of Normalized Sensory Nerve Axon Density: Data from different KIC mice depict the sympathetic axonal density within various ROI in the pancreas at 6.5 weeks and 8 weeks of age, categorized as ASYMP, ADM, PanIN, or PDAC. The quantification involved 6 ASYMP, 6 ADM, 6 PanIN, and 5 PDAC ROIs (Guillot et al., 2022). The dotted blue points correspond to the data. The green curve represents the mean density in each ROI, with identified stages on the x-axis extrapolated from the chronological knowledge of each ROI type.</figcaption></figure>
- The parameter $s_r$ determines the rate at which proliferation takes place. For small values of $s_r$, then the proliferation increases in rate in precancerous lesions, whereas for big values of $s_r$, the tumor growth increases in rate only in late stages of precancerous lesions. We choose $s_r \in [0.2, 10]$.

<figure><img src="image_16.png" /><figcaption>Figure 7: Illustration of the behaviour of the basal growth rate of the pre-tumor and tumor cells.</figcaption></figure>

* The growth rate of the pre-tumor and tumor cells can be regulated by the presence of axons through the term
  $-\mu_1 A_1(t) + \mu_2 A_2(t)$.

In that, we consider that the sympathetic axons either slow down or accelerate tumor growth whereas the sensory axons only accelerate tumor growth. The parameters $\mu_1$ and $\mu_2$ modulate the impact that the axons have on the proliferation. We consider that $\mu_1 \in [-1, 1]$ and $\mu_2 \in [10^{-4}, 1]$.

* It's important to point out that even if cells proliferate, the total quantity of cells is limited and there is no explosion in finite time. In other words, we require $1 - \frac{N(t)}{\tau_C} - \mu_1 A_1(t) + \mu_2 A_2(t) < \infty$ with $\tau_C$ being the carrying capacity of the amount of cells in the absence of axons. Thus, using the fact that $A_i \le 1$ for $i = 1, 2$ we have,

$$N(t) \le \tau_C(1 - \mu_1 A_1(t) + \mu_2 A_2(t)) \le \tau_C(1 + \mu_2 + \mu_1^{-}), \\ \le C(\tau_C),$$

with $C(\tau_C) = \tau_C(1 + \mu_2 + \mu_1^{-})$ where $x^{-} = \frac{1}{2}(|x| - x)$. Since $N(0) = 100$, we then require $N(0) \le \tau_C$. We choose $\tau_C \in [120, 300]$.

**Sympathetic axons growth.** We introduced an Allee effect in the growth term of the sympathetic axons (similar to what is seen in [Lolas et al., 2016]). We recall that the dynamics of $A_1$ is given by the following

$$r_{A_1} A_1(t) \left( \frac{A_1(t)}{\theta \left( \frac{N(t)}{N(0)} \right)} - 1 \right) (1 - A_1(t)).$$

* The parameter $r_{A_1}$ is the growth rate of the autonomic axons. The bigger the value of $r_{A_1}$ is, the steeper the increase or decrease of $A_1$ is. We consider $r_{A_1} \in [10^{-4}, 5]$.

* In the light of what is done in [Lolas et al., 2016], we also propose to modulate the threshold $\theta$ with the total amount of cells:

$$\theta \left( \frac{N(t)}{N(0)} \right) = \frac{A_1^{eq}}{2} + \frac{1}{2} \tanh \left( s_{\theta} \left( \frac{N(t)}{N(0)} - 1 - 0.1 \right) \right) + \frac{1}{2} \quad (11)$$

where $A_1^{eq}$ is the density of the sympathetic axons at the healthy state and $s_{\theta}$ is the steepness of the decrease speed of $A_1$. The value of $A_1^{eq}$ is estimated from the data (cf. Table 1). In order to ensure a biological meaning, the following has to hold: for small time $t$, $\theta(\frac{N(t)}{N(0)}) \approx \theta(1)$, and we want $\theta(1) < A_1(0)$ with $A_1(0) \approx A_1^{eq}$ to observe initially an increase in the autonomic axon density. Thus, we want

$$s_{\theta} \ge \frac{\tanh^{-1}(1 - A_1^{eq})}{10^{-1}}.$$

It leads us to consider $s_{\theta} \in [13, 20]$. For large values of the ratio $\frac{N(t)}{N(0)}$, we have $\theta(\frac{N(t)}{N(0)}) \approx 1 + \frac{A_1^{eq}}{2} > 1$ leading to a decrease in the autonomic axon density.

**Sensory axons growth.** The sensory axons follow a logistic growth regulated by the presence of PDAC

$$r_{A_2} \left( \frac{N_c(t)}{N(t)} \right) A_2(1 - A_2)$$

with $r_{A_2}(s) = \bar{r}_{A_2} \tanh(s_{A_2}s)$. The parameter $\bar{r}_{A_2}$ is the maximum growth rate of the sensory axons. We consider $\bar{r}_{A_2} \in [10^{-4}, 5]$. The parameter $s_{A_2}$ modulates the steepness of the increase of the growth of $A_2$. Here, we consider $s_{A_2} \in [10^{-4}, 10]$.
<figure><img src="image_20.png" /><figcaption>Fig. 11. Whisker Plot of Normalized Cancer Cell Concentration: Sympathectomy was conducted through the injection of neurotoxin 6-hydroxydopamine (6-OHDA) vehicle solution between 3.5 and 4 weeks of the KIC age. We utilized data published in Guillot et al. (2022), where the concentration of cancer cells in 6 KIC mice was quantified at 6.5 weeks. The mice were distributed as follows: 3 mice injected with ascorbic acid (AA) and 3 mice injected with the neurotoxin 6-hydroxydopamine (OHDA). For each mouse, we obtained 4 biopsies, and for each biopsy, we collected the proportion of PDAC. Therefore, we consider a total of 6 × 4 = 24 quantifications of the proportion of PDAC, with 3 × 4 = 12 quantifications from mice injected with AA and 12 quantifications from mice injected with OHDA. The dotted blue and red points correspond to the data for control mice and denervated mice respectively.</figcaption></figure>

The first coordinate of our criterion is defined as the mean squared error (MSE) given by

$$G_1(\theta) = \frac{2}{n} \sum_{i=1}^{n} \frac{1}{|I_{t_{\text{PDAC}}}|} \int_{I_{t_{\text{PDAC}}}} \left( p_i - p \left( t | \theta, x_i \right) \right)^2 dt \\ + \frac{1}{m} \sum_{j=n+1}^{n+m} \frac{1}{|I_j|} \int_{I_j} \left( a_j - A \left( t | \theta, z_j \right) \right)^2 dt. \quad (19)$$

In the first part of the MSE defined in (19), we have integrated the square difference over a time interval of size 6 because the exact time (on our time axis) at which the cell distributions were observed is not known

* $p(t|\theta, u_i) := N_c(t|\theta, u_i) / N(t|\theta, u_i)$ is the proportion of proliferative cells given by the model at time $t$, $N_c(t|\theta, u_i)$ being the amount of proliferative cells and $N(t|\theta, u_i)$ being the total amount of cells at time $t$ given by the model of Section 3.2.

* $I_{t_{\text{PDAC}}} = [t_{\text{PDAC}} - 3; t_{\text{PDAC}} + 3]$ is the time range at which the cell distributions were observed.

In the second part of the MSE defined in (19), we have integrated the square difference over a time interval of size 6 because the exact time (on our time axis) at which the axon densities were observed is not known, but we know the kind of tissues we sampled. More precisely,

* $I_j \in \{[t_{\text{ADM}} - 3; t_{\text{ADM}} + 3], [t_{\text{PANIN}} - 3; t_{\text{PANIN}} + 3], [t_{\text{PDAC}} - 3; t_{\text{PDAC}} + 3], [t_{\text{PDAC}}^{\text{advanced}} - 3; t_{\text{PDAC}}^{\text{advanced}} + 3]\}$ is the time range corresponding to the observation of the jth biological sample, i.e the sample coming from either the ADM, PanIN, PDAC or mature PDAC lobules.

* $A(t|\theta, z) := A_1(t|\theta)$ or $A_2(t|\theta)$ is the axon density given by the model at time $t$ of type $z$ and parameter set $\theta$.

The second component of our criterion is a penalization term to capture the correct behavior of the model according to the biological knowledge we have on the chronology of appearances of PDAC cells(see Aguirre et al., 2003) in the control group (AA). We do not expect many PDAC cells at time $t < t_{\text{PDAC}}^{\text{early}}$ where $t_{\text{PDAC}}^{\text{early}}$ is the first time

of appearance of PDAC lobules. Thus, the second component of our criterion is defined as

$$G_2(\theta) = \frac{1}{t_{\text{PDAC}}^{\text{early}}} \int_{0}^{t_{\text{PDAC}}^{\text{early}}} p(t|\theta, \text{AA}) dt, \quad (20)$$

where $p(t|\theta)$ is the proportion of proliferative cells at time $t$ given by the model. The quantity $G_2(\theta)$ is the mean proportion of proliferative cells in the control group (AA) over the time range $[0, t_{\text{PDAC}}^{\text{early}}]$ and it should be low to respect the biological knowledge we want to introduce in the calibration process.

### 4.3. Seeking for parameter sets with low criterion values

The range of parameters that makes sense in the mathematical model is the hypercube $\mathcal{H}$ given in Table 1. Without constraining $\theta \in \mathcal{H}$ with the above criterion $G(\theta)$, the dynamics of the model can change strongly based on the value of $\theta$. Many roads are available to restrain this range and obtain relevant parameter sets based on the 2-dimensional criterion $G(\theta) = (G_1(\theta), G_2(\theta))$. We could have tried to optimize a 1-dimensional criterion such as $\tilde{G}_{\lambda}(\theta) = G_1(\theta) + \lambda G_2(\theta)$, where $\lambda$ is a tuning parameter that defines the trade off between the data and the chronological information. Yet tuning $\lambda$ is difficult. Moreover, because of non-convexity, many local minima may exist, exhibiting different dynamics of the model. Instead of minimizing a $\tilde{G}_{\lambda}(\theta)$, we use a multi-stage algorithm that selects first the best sets $\theta$ of parameters according to $G_1(\theta)$ and then refines the selection according to $G_2(\theta)$.

The results given by the multi-stage algorithm given in Appendix B is a (relative large) collection of parameter sets $\theta$ that are consistent with the cell data and the chronological knowledge, i.e. that have low values of both $G_1(\theta)$ and $G_2(\theta)$. We start with a massive quasi-Monte Carlo sampling the hypercube $\mathcal{H}$ with a uniform distribution that acts here as a non-informative prior. The massive collection is then filtered according to both components of $G(\theta)$ and re-sample the hypercube $\mathcal{H}$ with an instrumental Gaussian distribution fitted on the filtered collection. The filtering step is then repeated to get the final collection of parameter sets.
## 3.2 Modeling denervation treatment

Testing the effects of denervation in preclinical models often requires the use of numerous animals, raising ethical concerns. Therefore, the use of a mathematical model to predict the effects of denervation becomes highly valuable, allowing for the replacement and reduction of the number of animals used. In the literature, the role of axons in cancer progression has been studied through chemical or surgical denervation. In our mathematical model, denervation is implemented by setting parameters related to the targeted axon to 0 after the time of the intervention.

Specifically, denervating sympathetic axons is equivalent to suppressing the influence of $A_1$ on cell transport and proliferation. This is achieved by setting the parameters as follows after the denervation time:

$$ \beta = 0, \mu_1 = 0. $$

Conversely, denervating sensory axons is equivalent to suppressing the influence of $A_2$ on cell transport and proliferation, after the denervation time, the parameters are set as:

$$ \delta = 0, \mu_2 = 0. $$

For denervating both axons, the parameters are set as follows after the denervation time:

$$ \beta = 0, \mu_1 = 0, \delta = 0, \mu_2 = 0. $$

We have to adapt the notations given above to take into account the above mathematical denervation when useful. The dynamical system defined in Section 2, in particular in Eq. (3) is parameterized by $\theta$ as explained in Section 3.1 and models the control treatment (AA). The total amount of cells and the total amount of cancer cells arising from this dynamical system is denoted by $N(t|\theta, AA)$ and $N_c(t|\theta, AA)$, respectively. The dynamical system that models the denervation treatment is obtained from the same dynamical system, yet the value of $\theta$ is changed after time of the intervention, where some components are set to 0 as explained above. The total amount of cells and the total amount of cancer cells arising from this second dynamical system is denoted by $N(t|\theta, OHDA)$ and $N_c(t|\theta, OHDA)$, respectively.

# 4 Calibration

In this section, we focus on the calibration of the model to the biological data at our disposal. It is rather complex in our case since

* a small number of biological data is provided, with variability in samples as described in Section 4.1,

* the model is far from being linear in the parameters and the predictions returned by the mathematical model cannot be written explicitly as a function of the parameters,

* we do not have at our disposal a likelihood function to explain the difference between the predictions of model and the data.

Moreover, the set of parameters is of rather high dimension ($d = 14$) and we need to avoid selecting the best set of parameters that would overfit the few observed data and would not reflect the biological phenomenon we are trying to capture.

To this aim, we propose in Section 4.2 a 2-dimension criterion to measure the quality of the calibration of the model to the data and to our biological knowledge explained in Section 4.1: (1) a mean squared error criterion based on the cell distribution and axon evolution, that measures the difference between model outcomes and data, and (2) a criterion based on the biological knowledge we have on the chronology of cell evolution. As described in Section 4.3, we managed to extract biologically relevant sets of parameters that are consistent with the cell data and the chronological knowledge by exploring the parameter space with a multi-stage quasi-Monte Carlo algorithm that starts from a flat, uniform prior over the parameter space which is the hypercube given by the range column of Table 1, and ends with the 0.4% best sets of parameters according to the criterion, that are presented in Section 4.4 and in Appendix B.2.

## 4.1 Biological data

The selection of experimental data for calibrating the mathematical model is crucial and should align with the characteristics of the relevant preclinical model. The calibration process encompasses various aspects of the animal model, incorporating chronological insights into cancer progression derived from the literature. Additionally, data on sympathetic nerve axons density, sensory nerve axons density, and cell density at different stages of cancer progression are considered during the calibration process.
<figure><img src="image_21.png" /><figcaption>Fig. 12. Schematic view of the step-wise algorithm used to select parameter sets. The collection $J_0$ of parameter sets is distributed uniformly over the hypercube $H$ of Table 1. The collection $J_1$ is obtained by filtering $J_0$ with $G_1(\theta) \le g_1$ and then $G_2(\theta) \le g_2$. The collection $J_2$ is obtained by re-sampling $H$ with an instrumental Gaussian distribution fitted on $J_1$. The collection $J_3$ is obtained by filtering $J_2$ with $G_1(\theta) \le g'_1$ and then $G_2(\theta) \le g'_2$. The thresholds $g_1$ and $g'_1$ are chosen as some quantiles of small order of the collection that is constrained. The final collection $J_3$ is a collection of parameter sets that are consistent with the cell data and the chronological knowledge, i.e. that have low values of both $G_1(\theta)$ and $G_2(\theta)$ (see Appendix B.2 for more details).</figcaption></figure>

## 4.4. Relevant sets of parameters

By taking some of the sets of parameters obtained in Section 4.3, as seen in Fig. 13, $A_1$ and $A_2$ have variability in their time of remodeling, as well as there is a difference in the time of arrival of cancer cells. We observe that despite the variability in dynamics, all parameter sets obtained fit the biological data. In the subsequent sections, we will utilize some of these parameter sets to conduct a detailed investigation of in silico denervation

# 5. Impact of the in silico denervation

The aim of this section is to investigate the effect of denervation on cancer progression using specific parameter sets selected from Section 4.4, chosen for their cost-effectiveness and alignment with biological data outlined in Table 2. Their values for the criteria of fitting the data (19) are 0.264, 0.2461, and 0.246 respectively. Their values for the criteria of fitting the chronological assumption (20) are 0.033, 0.002, and 0.014 respectively.

The parameter configurations 2 and 3 outlined in Table 2 are selected so that sympathetic axons induce an inhibitory effect on tumor growth ($\mu_1 > 0$). However, they are distinguished by varying durations for axon remodeling, reflecting patterns observed in the literature (Guillot et al., 2022) and accommodating the variability inherent in biological data. Furthermore, in light of findings from a literature source (Renz et al., 2018a) suggesting a positive impact of $A_1$ on tumor growth, we also choose parameter set 1 with $\mu_1 < 0$ to investigate this case.

In what follows, we analyze the effects of varying denervation timings on either sympathetic or sensory axons, or both, in order to elucidate their impacts on tumorigenesis.

## 5.1. Indicator of the invasive potential

To study the effect of denervation and the role of axons in pancreatic cancer progression, we investigate the effect of three types of denervation:

* the denervation of sympathetic axons ($A_1$),

* the denervation of sensory axons ($A_2$),

* the denervation of both types of axons ($A_1$ and $A_2$).

A quantitative measure of the effect of denervation for any of the above three types of denervation can be the difference at a given time between the amount of cancer cells in the control case (with both axons regulating the dynamics) and in the denervated case. Although this function allows us to measure the effect of denervation, it does not allow us to take into account all the specificities of temporal dynamics. Indeed, the earlier a cancer cell appears, the higher the probability of accumulating genetic alterations that enhance aggressiveness and metastasis. So we construct the following function to measure the impact of denervation:

$$I(T) = \frac{1}{T} \left( \int_{0}^{T} \frac{N_c(t|\theta, \text{OHDA})}{N(t|\theta, \text{OHDA})} dt - \int_{0}^{T} \frac{N_c(t|\theta, \text{AA})}{N(t|\theta, \text{AA})} dt \right) \quad (21)$$

with the notations of Section 3.2. This indicator gives information about the invasive potential at time $T$ because it takes into account the history of existence of cancer cells over time.

This indicator is positive for an overall pro-tumoral effect of denervation ($I(T) > 0$), and conversely it is negative for an overall anti-tumoral effect of denervation ($I(T) < 0$).

## 5.2. In silico denervation at defined time points

In this section, we apply the three types of denervations stated in Section 5.1. These denervations are implemented at two specific times:

* at 28 days, which corresponds to the time between the observation of PanIN and the onset of PDAC for the KIC model (cf. Section 4.1),

* at 50 days, which corresponds to the time between the onset of PDAC and the observation of advanced PDAC (cf. Section 4.1).

The dynamics of the model subject to these denervations are illustrated in Fig. 14.

We denote by $T$ the time of observation, then the effects of the denervations are quantified by the indicator $I(T)$ given by (21) and the results are summarized in Table 3. The invasive potential can be quantitatively compared with the various results in the literature on the effect of axons on the initiation and progression of PDAC:

* in Saloman et al. (2016), the authors show that ablating sensory axons prior to the onset of the pathology (referred to as early denervation) inhibits tumorigenesis of the pancreatic cancer. This result is translated by $I(T) < 0$.

* In Guillot et al. (2022), the authors show that ablating sympathetic axons prior to the onset of the pathology (referred to as early denervation) leads to an acceleration of tumor growth and metastasis. This result is translated by $I(T) > 0$ in the case of an early denervation of $A_1$.

* Finally, in Renz et al. (2018a) and Guillot et al. (2022), the authors conduct surgical denervation for both $A_1$ and $A_2$ after and before invasive tumor formation, categorized as early and late denervation, respectively. This procedure results in the ablation of mixed sympathetic and sensory axons.

- In Renz et al. (2018a) they show that removing the axons inhibits tumor growth. Hence, this result is that $I(T) < 0$ in the case of late denervation of $A_1$ and $A_2$.

- However, in Guillot et al. (2022), early denervation of the splanchnic nerve leads to a pro-tumoral effect with observations of metastasis and a smaller survival time of mice. This result is translated by $I(T) > 0$ in the case of early denervation of $A_1$ and $A_2$.
**Mice model and data acquisition methods.** In this study, we use the Kras<sup>LSL-G12D/+</sup>; Cdkn2a(Ink4a/Arf)<sup>lox/lox</sup>; Pdx1-Cre (KIC) mouse model of pancreatic cancer. The over-expression of a mutated form of the oncogene KRAS in pancreatic cells induces transdifferentiation, where healthy acinar cells (the functional unit of the exocrine pancreas) transform into ductal cell-like cells, leading to the formation of acinar to ductal metaplasia (ADM). Subsequently, ADM can progress to form premalignant pancreatic intraepithelial neoplasia (PanIN). The deletion of tumor-suppressor genes Ink4A and Arf further accelerates tumor development, resulting in PDAC itself. The KIC mice model successfully replicates the stepwise progression observed in human pancreatic cancer (cf. [Aguirre et al., 2003], and Figure 8). Several mouse pancreases at the time points of 6.5 weeks and 8 weeks are utilized in the experiments to accommodate biological heterogeneity. The technique used to quantify nerve axons density and cell density within lobules of the pancreatic tissue is detailed in the source data file of [Guillot et al., 2022]. In short, data comes from the quantification of 3D cleared tissue images obtained by Light Sheet Fluorescence Microscopy. Tissue staining was achieved through immunostaining and tissue clearing using the iDISCO+ protocol. Tyrosine hydroxylase antibody staining (TH), an enzyme involved in the biosynthesis of norepinephrine (one of the main neurotransmitters of the sympathetic nervous system), and Calcitonin gene-related peptide (CGRP) antibody were used to visualize sympathetic and sensory innervation respectively. Regions of Interest (ROI), namely asymptomatic (ASYM), ADM, PanIN, and PDAC, were segmented based on the autofluorescence signal of the tissue, and their volumes were measured. Axonal networks were manually reconstructed to measure the total axon length in each ROI. The unitless axon density was calculated as follows: ((axons length sum × 1000)(µm)/(volume of each Regions of Interest (ROI))<sup>1/3</sup> (µm)). Densities obtained through experimentation were normalized by first dividing by the cubic root of the volume in which the density was measured and then by an affine transformation to obtain proportions. This normalization facilitates the comparison of quantities both within the experimental data and with the outputs provided by the mathematical model (cf. Figures 9 and 10).

The method used to describe the proportion of tissue categorized by their phenotypes in histological sections through the pancreas of a 6.5-week-old mouse can be found in the source data file of [Guillot et al., 2022] (cf. Figure 11).

**Knowledge on the chronological process.** Based on the previous characterization of the KIC model (seen in [Aguirre et al., 2003] and in our personal observation), the time of observation of ADM lobules denoted $t_{ADM}$ is around 17 days, the first appearance of relatively advanced-staged PanIN lobules denoted $t_{PANIN}$ is observed around 24.5 days, and the first appearance of PDAC lobules denoted $t_{PDAC}^{early}$ is observed around 35 days. The median survival of a KIC mouse model is 9 weeks and so, at 6.5 weeks ($= t_{PDAC}$) the mice has already developed a tumor, but at 8 weeks ($= t_{PDAC}^{advanced}$), the cancer is really aggressive. All these chronological assumptions are summarized in Figure 8. Additionally, we include a time interval of plus or minus three days around the above temporal values to account for variability in observations.

<figure><img src="image_17.png" /></figure>

Figure 8: Diagram outlining the chronological assumptions regarding cancer progression in pancreas of KIC mice, drawing from [Aguirre et al., 2003] and our firsthand observations. The indicated times represent the estimated average time at which each new stage is expected to be observed for the first time. The bold red segments indicate a time interval of plus or minus three days around the mentioned temporal values, encompassing variations in chronological estimates.

**Sympathetic and Sensory nerve axons quantification.** The normalized densities of sympathetic and sensory nerve axons are shown in Figure 9 and Figure 10, respectively, quantified at 6.5 weeks of the KIC age. For sensory axons, density measurements were also performed at 8 weeks for the PDAC lesions. Also, each sample of density measurements from the data does not only characterize the axon densities but also the tissue stages from where it was observed (ROI). This additional information may not be fully exploited by the model because information about space is not taken into account in the construction of the model. One way to use all information of the data is to extrapolate a dynamic based on the ROI type of the data origin. Since we know that $t_{ASYM}$, $t_{ADM}$, $t_{PANIN}$, and $t_{PDAC}$ are the times of observation of each ROI type during cancer progression, we then consider that the quantification of the normalized densities of the axons in each ROI type corresponds to those times respectively. Those timings are stated on the x-axis of the Figure 9 and Figure 10. Hence, the dynamics of the axons are illustrated by the green curves in Figures 9 and 10. It allows to combine qualitative information (chronological knowledge) to quantitative information and ensures a more precise calibration of the mathematical model.
<figure><img src="image_22.png" /><figcaption>**Fig. 13.** Dynamics of the model using different sets of parameters from the relevant sets of parameters obtained. (Top) Evolution of the sympathetic axons ($A_1$). (Middle) Evolution of the sensory axons ($A_2$). (Bottom) Evolution of the proportion of cancer cells ($N_c/N$). The blue dots in each sub-figure correspond to the biological data provided for the sympathetic axons, sensory axons, and the frequency of cancer cells respectively. The red dots in the bottom frame correspond to the frequency of cancer cells observed at 6.5 weeks after an early denervation at 28 days with the use of OHDA, and the dashed red curves correspond to the evolution of the concentration of cancer cells after early denervation.</figcaption></figure>

<table><caption>**Table 2**

Different sets of parameters.</caption><thead><tr><th></th><th>π<sub>0</sub></th><th>β</th><th>δ</th><th>γ<sub>r</sub></th><th>s<sub>r</sub></th><th>τ<sub>C</sub></th><th>μ<sub>1</sub></th><th>μ<sub>2</sub></th><th>r<sub>A<sub>1</sub></sub></th><th>r&#773;<sub>A<sub>2</sub></sub></th><th>x<sub>1,π</sub></th><th>ε<sub>1,π</sub></th><th>s<sub>θ</sub></th><th>s<sub>A<sub>1</sub></sub></th></tr></thead><tbody><tr><td>Set 1</td><td>2.005</td><td>0.73</td><td>0.398</td><td>1.50</td><td>2.68</td><td>172.295</td><td>-0.176</td><td>0.214</td><td>0.055</td><td>0.241</td><td>32.97</td><td>5.357</td><td>15.131</td><td>4.151</td></tr><tr><td>Set 2</td><td>4.589</td><td>0.504</td><td>0.829</td><td>1.160</td><td>2.775</td><td>177.807</td><td>0.609</td><td>0.139</td><td>0.077</td><td>0.928</td><td>34.56</td><td>6.555</td><td>14.733</td><td>1.105</td></tr><tr><td>Set 3</td><td>1.795</td><td>0.535</td><td>0.398</td><td>4.537</td><td>6.097</td><td>150.576</td><td>0.176</td><td>0.678</td><td>0.032</td><td>0.29</td><td>30</td><td>4.385</td><td>17.609</td><td>6.49</td></tr></tbody></table>

<table><caption>**Table 3**

Values in percentage of the invasive potential $I(70)$ (cf. (21)) for the sets of parameters of Table 2. The columns indicate the different denervations, sympathetic axons ($A_1$), sensory axons ($A_2$) or both types of axons ($A_1$ and $A_2$), at early stage (28 days) or late stage (50 days).</caption><thead><tr><th></th><th>A<sub>1</sub> at 28d</th><th>A<sub>1</sub> at 50d</th><th>A<sub>2</sub> at 28d</th><th>A<sub>2</sub> at 50d</th><th>A<sub>1</sub> and A<sub>2</sub> at 28d</th><th>A<sub>1</sub> and A<sub>2</sub> at 50d</th></tr></thead><tbody><tr><td>Set 1</td><td>-2.311</td><td>-0.56</td><td>-1.697</td><td>-1.58</td><td>-3.961</td><td>-2.249</td></tr><tr><td>Set 2</td><td>5.893</td><td>0.876</td><td>-2.163</td><td>-1.85</td><td>3.435</td><td>-0.788</td></tr><tr><td>Set 3</td><td>2.493</td><td>0.441</td><td>-6.94</td><td>-6.366</td><td>-4.482</td><td>-5.47</td></tr></tbody></table>

The only nuance to note is that the time at which the effect of denervations is observed in the various articles is not so easily translated into a quantifiable datum $T$ in days. For instance, the mice models used may differ in terms of chronological progression of the disease and the different experimental techniques used may also affect the temporal component of the data in their own way. To make the most of this variability, the final time of observation $T$ is set at $T = 70$ days as the rounded upper bound of the median survival of mice (63 days in Guillot et al., 2022) which allows us to observe the effect of the denervation over the entire time range.

### 5.2.1. Case when $\mu_1 < 0$

Based on Table 3, the invasive potential of parameters set 1 is always negative for the three types of denervations. Indeed, only an anti-tumoral result is observed for the early and late denervation as seen in Fig. 14(b), which do not recapitulate the results shown in Guillot et al. (2022). Therefore, this set can be considered as a false positive estimation due to the high variability of the data and the ill-posedness of the calibration problem. Thus, in what follows, we focus our study on the parameters sets 2 and 3 of Table 2 when $\mu_1 > 0$.

### 5.2.2. Case when $\mu_1 > 0$

In this section, we analyze the impact of axons on tumorigenesis using parameters sets 2 and 3 from Table 2. Our model posits that sympathetic axons exert an anti-tumoral effect. Therefore, denervation of $A_1$ leads to an increased proportion of cancerous cells revealed by a positive value for invasive potential, whether denervation occurs early or late (see Fig. 14 and Table 3). Conversely, sensory axons are implicated in promoting tumorigenesis. Thus, denervation of $A_2$ leads to a decrease in proportion of cancerous cells and a negative value for invasive potential (see Fig. 14 and Table 3).

**Effect of the double denervation of sympathetic and sensory axons.** The impact of denervating sympathetic and sensory axons proves to be complex and subject to diverse outcomes across the two distinct parameter sets, reflecting biological variability.

**Outcome 1.** For both set of parameters, early denervation of sympathetic and sensory axons result in a shift from an initial pro-tumoral effect of denervation (increase proportion of cancer cell compare to control condition) to a late anti-tumoral effect of denervation (decrease proportion of cancer cell compare to control condition). This “shift
<figure><img src="image_18.png" /><figcaption>Figure 9: Whisker Plot of Normalized Sympathetic Nerve Axon Density: Data from different KIC mice depict the sympathetic axonal density within various ROI in the pancreas at 6.5 weeks of age, categorized as ASYMP, ADM, PanIN, or PDAC. The quantification involved 4 ASYMP, 5 ADM, 6 PanIN, and 7 PDAC ROIs [Guillot et al., 2022]. The dotted blue points correspond to the data. The green curve represents the mean density in each ROI, with identified stages on the x-axis extrapolated from the chronological knowledge of each ROI type.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 10: Whisker Plot of Normalized Sensory Nerve Axon Density: Data from different KIC mice depict the sympathetic axonal density within various ROI in the pancreas at 6.5 weeks and 8 weeks of age, categorized as ASYMP, ADM, PanIN, or PDAC. The quantification involved 6 ASYMP, 6 ADM, 6 PanIN, and 5 PDAC ROIs [Guillot et al., 2022]. The dotted blue points correspond to the data. The green curve represents the mean density in each ROI, with identified stages on the x-axis extrapolated from the chronological knowledge of each ROI type.</figcaption></figure>

**PDAC Cell Concentration.** Figure 11 illustrates the data extracted for the normalized cancer cell concentration at 6.5 weeks. We set the time of observation for the data and aim to identify model dynamics that closely match these data points at 6.5 weeks. Through observation, we note that the proportion of PDAC cells is higher in OHDA samples than in AA samples.

## 4.2 A 2-dimensional criterion

The first dimension of our criterion is a mean squared error criterion that measures the difference between model outcomes and data. We consider the following data:
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 14.** Model dynamics subject to denervation. (Left) Evolution of the axons. (Right) Evolution of the proportion of cancerous cells $\frac{N_c(t)}{N(t)}$. The different colors show the effect of different denervations. The first row (resp. second row and third row) corresponds to the dynamics for the parameters set 1 (resp. 2 and 3) in Table 2.</figcaption></figure>

effect” is highlighted by the crossing of the green denervated and the black control curves in Fig. 14.

**Outcome 2.** For both set of parameters, late denervation of both axons types inhibits tumorigenesis. In Fig. 14, the green denervated curves are below the control ones meaning that pathological cancerous cells are present at a lower level compared to the control scenario.

**Outcome 3.** However significant difference emerges between the two sets. For early denervation, invasive potential is positive for parameter set 2 and negative for parameter set 3. This suggests a higher net

production of cancerous cells over the entire time span for set 2 compared to set 3 (higher area between the green and black curves before the “shift effect” and a lower area after the curves crossing, see Fig. 14). This difference between the two sets reflects the dominance of a specific axon type (sympathetic for set 2 and sensory for set 3) in pancreatic cancer regulation.

The role of axons in cancer progression can be effectively established using the mathematical model and in silico denervation experiments. However, the relevance of in silico results is partly correlated with the calibration of the model. This calibration becomes more
<figure><img src="image_20.png" /><figcaption>Figure 11: Whisker Plot of Normalized Cancer Cell Concentration: Sympathectomy was conducted through the injection of neurotoxin 6-hydroxydopamine (6-OHDA) vehicle solution between 3.5 and 4 weeks of the KIC age. We utilized data published in [Guillot et al., 2022], where the concentration of cancer cells in 6 KIC mice was quantified at 6.5 weeks. The mice were distributed as follows: 3 mice injected with ascorbic acid (AA) and 3 mice injected with the neurotoxin 6-hydroxydopamine (OHDA). For each mouse, we obtained 4 biopsies, and for each biopsy, we collected the proportion of PDAC. Therefore, we consider a total of 6 × 4 = 24 quantifications of the proportion of PDAC, with 3 × 4 = 12 quantifications from mice injected with AA and 12 quantifications from mice injected with OHDA. The dotted blue and red points correspond to the data for control mice and denervated mice respectively.</figcaption></figure>

* n biological samples of cell distributions $p_i \in (0; 1)$ under treatment $u_i \in \{AA, OHDA\} (i = 1, \dots, n)$ that should be compared to the proportion of proliferating cells $p(t|\theta, u)$ given by the model, on average over time $t \in (0, T)$

* m biological samples of axon density $a_i \in (0; 1)$ in AA ($i = n + 1, \dots, n + m$) that should be compared to the axon density evolution $A(t|\theta)$ given by the model, on average over the time range $t \in I_i$ corresponding to the kind of tissues of i-th biological sample (ADM, PanIN, PDAC or mature PDAC). Moreover, $z_i \in \{A_1, A_2\}$ indicates whether the observed density $a_i$ is a density of sympathetic ($A_1$) or sensory ($A_2$) axons.

The first coordinate of our criterion is defined as the mean squared error (MSE) given by

$$G_1(\theta) = \frac{2}{n} \sum_{i=1}^{n} \frac{1}{|I_{t_{\text{PDAC}}}|} \int_{I_{t_{\text{PDAC}}}} \left( p_i - p(t|\theta, x_i) \right)^2 dt + \frac{1}{m} \sum_{j=n+1}^{n+m} \frac{1}{|I_j|} \int_{I_j} \left( a_j - A(t|\theta, z_j) \right)^2 dt. \quad (12)$$

In the first part of the MSE defined in (12), we have integrated the square difference over a time interval of size 6 because the exact time (on our time axis) at which the cell distributions were observed is not known

* $p(t|\theta, u_i) := N_c(t|\theta, u_i) / N(t|\theta, u_i)$ is the proportion of proliferative cells given by the model at time $t$, $N_c(t|\theta, u_i)$ being the amount of proliferative cells and $N(t|\theta, u_i)$ being the total amount of cells at time $t$ given by the model of Section 3.2.

* $I_{t_{\text{PDAC}}} = [t_{\text{PDAC}} - 3; t_{\text{PDAC}} + 3]$ is the time range at which the cell distributions were observed.

In the second part of the MSE defined in (12), we have integrated the square difference over a time interval of size 6 because the exact time (on our time axis) at which the axon densities were observed is not known, but we know the kind of tissues we sampled. More precisely,

* $I_j \in \{[t_{\text{ADM}} - 3; t_{\text{ADM}} + 3], [t_{\text{PANIN}} - 3; t_{\text{PANIN}} + 3], [t_{\text{PDAC}} - 3; t_{\text{PDAC}} + 3], [t_{\text{PDAC}}^{\text{advanced}} - 3; t_{\text{PDAC}}^{\text{advanced}} + 3]\}$ is the time range corresponding to the observation of the j-th biological sample, i.e the sample coming from either the ADM, PanIN, PDAC or mature PDAC lobules.

* $A(t|\theta, z) := A_1(t|\theta)$ or $A_2(t|\theta)$ is the axon density given by the model at time $t$ of type z and parameter set $\theta$.

The second component of our criterion is a penalization term to capture the correct behavior of the model according to the biological knowledge we have on the chronology of appearances of PDAC cells(see [Aguirre et al., 2003])
<figure><img src="image_29.png" /><figcaption>Fig. 15. Evolution of the invasive potential with respect to varying denervation times for the parameters set 2 from Table 2. The times of denervations of $A_1$ (resp. $A_2$) are indicated on the y-axis (resp. the x-axis). The observation time $s_i$ of the invasive potential is indicated in days above each heatmap.</figcaption></figure>

complex when the data is scarce and highly variable. Consequently, the parameter sets obtained during calibration exhibit different dynamics and denervation effects. Moreover, it is important to note that both the timing of denervation and the duration of observation post-denervation are critical factors in understanding the denervation's effects.

### 5.3. In silico denervation for varying times

Using the mathematical model, we investigate the temporal dynamics by denervating both sympathetic and sensory axons at various time points. Subsequently, we construct heatmaps at different observation times to have an evolving illustration of the invasive potential.

#### 5.3.1. Heatmap construction and interpretation

Studying experimentally the effect of different types of denervation can be time consuming and expensive. However, the mathematical model allows to perform a large number of in silico denervations and hence gives insights on the role of the axons in tumorigenesis. In this Section, using the parameters sets 2 and 3 of Table 2, we perform denervations varying the time of denervation of both axons independently and we quantify their pro- or anti-tumoral effect at different times during the tumorigenesis process. The results of the in silico denervations are illustrated by the two sequences of heatmaps of Figs. 15 and 16.

First, we introduce the finite sequence $(s_i) \in \{5, 10, ..., 70\}$ which corresponds to the observation times for the effect of denervations. Then, we denote by $t_{A_1}$ (resp. $t_{A_2}$) the time of denervation of the sympathetic axons $A_1$ (resp. the sensory axons $A_2$). Each heatmap of Fig. 15 or 16 corresponds to a grid where the y-axis corresponds to the different values taken by $t_{A_1}$ and the x-axis to those taken by $t_{A_2}$. The following holds

$$t_{A_j} \in [0, s_i], \quad j = 1, 2,$$

since it makes no sense to look at the effect of denervation before denervation has taken place.

The cell from the heatmap $s_i$ located at $(t_{A_1}, t_{A_2})$ corresponds to the invasive potential at time $s_i$ subject to the denervation of sympathetic axons at time $t_{A_1}$ and subject to the denervation of sensory axons at time $t_{A_2}$ (cf. (21) and Section 5.1). We keep track of the denervation

time in the notation of the index evaluated at time $s_i$ as follows: $I(s_i; t_{A_1}, t_{A_2})$. Hence, two kind of results are observed:

* if $I(s_i; t_{A_1}, t_{A_2}) > 0$, then the cell is colored in blue. Hence, the aggressiveness of the tumor is encoded by this color. The darker the blue, the stronger the pro-tumoral effect of the denervation.

* If $I(s_i; t_{A_1}, t_{A_2}) < 0$, then the cell is colored in brown. Hence, the possibility of tumor remission through denervation is encoded by this color. The darker the brown, the stronger the anti-tumoral effect of the denervation.

Moreover, the sequence of heatmaps illustrates the dynamical evolution of the invasive potential. For instance, the red, blue, green, and black stars on the heatmaps of Fig. 15 (resp. Fig. 16) correspond to the invasive potential of the early denervated and control curves of the parameters set 2 (resp. set 3) of Table 2 illustrated by the red, blue, green and black curves of Fig. 14(d) (resp. Fig. 14(f)).

#### 5.3.2. Impact of the time-varying denervations on the evolution of the tumorigenesis

In the following, a more qualitative approach of model validation is proposed. A detailed description of the interesting elements of the model dynamics subject to time-varying denervations is given. We discuss the possibility of establishing the best time and strategy of denervation.

Different times of remodeling of sympathetic axons explain the different apparition times of the pro-tumoral effect of denervations. The pro-tumoral effect of denervation starts to be observed at $s_i = 25$ days for parameters set 2 (cf. Fig. 15) and at $s_i = 40$ days for parameters set 3 (cf. Fig. 16) due to the darkening of the blue color in the lower halves of the heatmaps. The difference in the time of observation $s_i$ can be explained due to the difference in the time of arrival of sympathetic axons and the difference in the impact they have on cancer growth. By referring to Table 2, the parameters $r_{A_1}$ and $\mu_1$ are larger for set 2 than set 3, explaining the larger and earlier effect that sympathetic axons denervation exhibits on tumor growth as seen in Fig. 15 compared to Fig. 16.

In the early stages of the tumorigenesis, the strong pro-tumoral effect of denervation is associated to the early denervation of sympathetic axons. For $s_i \le 50$, the denervations associated with
in the control group (AA). We do not expect many PDAC cells at time $t < t_{\text{PDAC}}^{\text{early}}$ where $t_{\text{PDAC}}^{\text{early}}$ is the first time of appearance of PDAC lobules. Thus, the second component of our criterion is defined as

$$G_2(\theta) = \frac{1}{t_{\text{PDAC}}^{\text{early}}} \int_{0}^{t_{\text{PDAC}}^{\text{early}}} p(t|\theta, \text{AA})dt, \qquad (13)$$

where $p(t|\theta)$ is the proportion of proliferative cells at time $t$ given by the model. The quantity $G_2(\theta)$ is the mean proportion of proliferative cells in the control group (AA) over the time range $[0, t_{\text{PDAC}}^{\text{early}}]$ and it should be low to respect the biological knowledge we want to introduce in the calibration process.

## 4.3 Seeking for parameter sets with low criterion values

The range of parameters that makes sense in the mathematical model is the hypercube $\mathcal{H}$ given in Table 1. Without constraining $\theta \in \mathcal{H}$ with the above criterion $G(\theta)$, the dynamics of the model can change strongly based on the value of $\theta$. Many roads are available to restrain this range and obtain relevant parameter sets based on the 2-dimensional criterion $G(\theta) = (G_1(\theta), G_2(\theta))$. We could have tried to optimize a 1-dimensional criterion such as $\bar{G}_{\lambda}(\theta) = G_1(\theta) + \lambda G_2(\theta)$, where $\lambda$ is a tuning parameter that defines the trade off between the data and the chronological information. Yet tuning $\lambda$ is difficult. Moreover, because of non-convexity, many local minima may exist, exhibiting different dynamics of the model. Instead of minimizing a $\bar{G}_{\lambda}(\theta)$, we use a multi-stage algorithm that selects first the best sets $\theta$ of parameters according to $G_1(\theta)$ and then refines the selection according to $G_2(\theta)$.

The results given by the multi-stage algorithm given in Appendix B is a (relative large) collection of parameter sets $\theta$ that are consistent with the cell data and the chronological knowledge, i.e. that have low values of both $G_1(\theta)$ and $G_2(\theta)$. We start with a massive quasi-Monte Carlo sampling the hypercube $\mathcal{H}$ with a uniform distribution that acts here as a non-informative prior. The massive collection is then filtered according to both components of $G(\theta)$ and re-sample the hypercube $\mathcal{H}$ with an instrumental Gaussian distribution fitted on the filtered collection. The filtering step is then repeated to get the final collection of parameter sets.

<figure><img src="image_21.png" /><figcaption>Figure 12: Schematic view of the step-wise algorithm used to select parameter sets. The collection $\mathcal{J}_0$ of parameter sets is distributed uniformly over the hypercube $\mathcal{H}$ of Table 1. The collection $\mathcal{J}_1$ is obtained by filtering $\mathcal{J}_0$ with $G_1(\theta) \le g_1$ and then $G_2(\theta) \le g_2$. The collection $\mathcal{J}_2$ is obtained by re-sampling $\mathcal{H}$ with an instrumental Gaussian distribution fitted on $\mathcal{J}_1$. The collection $\mathcal{J}_3$ is obtained by filtering $\mathcal{J}_2$ with $G_1(\theta) \le g'_1$ and then $G_2(\theta) \le g'_2$. The thresholds $g'_i$ and $g'_i$ are chosen as some quantiles of small order of the collection that is constrained. The final collection $\mathcal{J}_3$ is a collection of parameter sets that are consistent with the cell data and the chronological knowledge, i.e. that have low values of both $G_1(\theta)$ and $G_2(\theta)$ (see Appendix B.2 for more details).</figcaption></figure>

## 4.4 Relevant sets of parameters

By taking some of the sets of parameters obtained in Section 4.3, as seen in Figure 13, $A_1$ and $A_2$ have variability in their time of remodeling, as well as there is a difference in the time of arrival of cancer cells. We observe that despite the variability in dynamics, all parameter sets obtained fit the biological data. In the subsequent sections, we will utilize some of these parameter sets to conduct a detailed investigation of in silico denervation
<figure><img src="image_30.png" /><figcaption>**Fig. 16.** Evolution of the invasive potential with respect to varying denervation times for the parameters set 3 from Table 2. The times of denervations of $A_1$ (resp. $A_2$) are indicated on the y-axis (resp. the x-axis). The observation time $s_i$ of the invasive potential is indicated in days above each heatmap.</figcaption></figure>

a strong pro-tumoral effect are localized in the lower halves of the heatmaps for both sets of parameters, corresponding to early denervation of sympathetic axons. The time of denervation of sensory axons has almost no effect in the lower halves of the heatmaps, as the blue region extends homogeneously along the x-axis. The fact that the denervations of sensory axons do not play a significant role in this time range can be explained due to the late remodeling of sensory axons and the delay between the denervation time of this type of axons and its impact on the dynamics of the model. Thus, in that time range, all types of denervations are mainly impacted by the denervation of the sympathetic axons.

In the late stages of the tumorigenesis, the impact of sensory axons' denervation becomes significant as the dynamics for both parameters sets undergo the shift effect described in the previous section. For $s_i \ge 55$, the sensory axons starts to display its inhibiting effect on cancer cells which can be seen through expanding brown areas or lightening blue areas on the upper halves of the heatmaps. It illustrates the antagonistic role of both types of axons: the sympathetic axons playing an anti-tumoral role on the early stages of tumorigenesis and the sensory axons playing a pro-tumoral role on the late stages of the tumorigenesis.

The sympathetic axons' (resp. sensory axons') impact on cancer growth dominates in Fig. 15 (resp. Fig. 16). On the one hand, the regulation of the sympathetic axons plays the most significant role on tumorigenesis for parameters set 2 since the denervation of this type of axons highlights a strong deviation of its dynamic from the control one (illustrated by the opacity level of the blue region in Fig. 15). On the other hand, the regulation of the sensory axons dominates for the parameters set 3 since the anti-tumoral effect of denervation becomes more and more significant as time passes. The brown area in Fig. 16 expands from the upper region at $s_i = 55$ days to three-quarter of the heatmap's area at $s_i = 70$ days (everywhere but the south-east corner corresponding to the latest denervation of sensory axons and the earliest denervation of sympathetic axons).

The effect of denervation takes time to be seen. Although an initial pro-tumoral effect of the denervation of both types of axons is always seen due to the dominant blue areas in the heatmaps for small/intermediate values of $s_i$, a later anti-tumoral effect of denervation is also observed when compared to the control curve after the first

remodeling of sensory axons. If the effect of sensory axons on cancer growth is strong enough, then the arrival of sensory axons may be sufficient for the overall anti-tumoral effect caused by it is denervation to compensate the previous pro-tumoral effect that has taken place. In Fig. 16, if we compare the heatmaps at $s_i = 50$ and $s_i = 60$ days, we see areas going from blue to brown. In that case, the final state of the disease (for large $s_i$) is a reduction of tumor, contrary to what can be observed initially for small $s_i$. In Fig. 15, the stronger pro-tumoral effect resulting from the denervation of $A_1$ takes longer to be compensated by the anti-tumoral effect of the denervation of $A_2$. However, it is possible to conjecture that the dynamics obtained of parameters set 2 at a later observation time ($s_i \gg 70$ days) will be similar to the ones displayed in the last frame Fig. 16 (parameters set 3 at observation time 70 days).

## 6. Conclusion

This study presents mathematical tools to model and simulate the joint effect of PNS axons (promoting and/or inhibiting cancer progression and proliferation) in pancreatic cancer tumorigenesis. It extends the previous model presented in Chauvet et al. (2023) by considering the cell phenotype as a continuous variable. This new mathematical formalism provides a more accurate description of tumor progression and associated neuroplastic changes. The mathematical model is then finely calibrated to the available data by measuring the goodness of fit between the model output and the biological knowledge with a two-dimensional criterion and by selecting the relevant calibrated parameter sets with a multi-stage quasi-Monte Carlo algorithm. In addition, a quantitative indicator of the balance between the two opposing pro- and anti-tumor effects of denervation is calculated numerically. This balance can be visualized over time for several parameter sets reflecting biological variability.

From a biological point of view, the mathematical model reconciles all the biological data found in the literature and provides an explanation for the opposite effects of surgical denervation performed at early and late stages of PDAC. Specifically, the mathematical model shows that when sympathetic axons increase tumor growth, the model does not replicate the biological data from the literature. This rules out the hypothesis of a functional switch of the sympathetic nervous system during tumorigenesis. However, when sympathetic axons have a consistently inhibitory effect on tumor growth, the model recapitulates all the
<figure><img src="image_22.png" /><figcaption>Figure 13: Dynamics of the model using different sets of parameters from the relevant sets of parameters obtained. (Top) Evolution of the sympathetic axons ($A_1$). (Middle) Evolution of the sensory axons ($A_2$). (Bottom) Evolution of the proportion of cancer cells ($N_c/N$). The blue dots in each sub-figure correspond to the biological data provided for the sympathetic axons, sensory axons, and the frequency of cancer cells respectively. The red dots in the bottom frame correspond to the frequency of cancer cells observed at 6.5 weeks after an early denervation at 28 days with the use of OHDA, and the dashed red curves correspond to the evolution of the concentration of cancer cells after early denervation.</figcaption></figure>

# 5 Impact of the in silico denervation

The aim of this section is to investigate the effect of denervation on cancer progression using specific parameter sets selected from Section 4.4, chosen for their cost-effectiveness and alignment with biological data outlined in Table 2. Their values for the criteria of fitting the data (12) are 0.264, 0.2461, and 0.246 respectively. Their values for the criteria of fitting the chronological assumption (13) are 0.033, 0.002, and 0.014 respectively.

The parameter configurations 2 and 3 outlined in Table 2 are selected so that sympathetic axons induce an inhibitory effect on tumor growth ($\mu_1 > 0$). However, they are distinguished by varying durations for axon remodeling, reflecting patterns observed in the literature [Guillot et al., 2022] and accommodating the variability inherent in biological data. Furthermore, in light of findings from a literature source [Renz et al., 2018a] suggesting a positive impact of $A_1$ on tumor growth, we also choose parameter set 1 with $\mu_1 < 0$ to investigate this case.

In what follows, we analyze the effects of varying denervation timings on either sympathetic or sensory axons, or both, in order to elucidate their impacts on tumorigenesis.

<table><caption>Table 2: Different sets of parameters</caption><thead><tr><th></th><th>π<sub>0</sub></th><th>β</th><th>δ</th><th>γ<sub>r</sub></th><th>s<sub>r</sub></th><th>τ<sub>C</sub></th><th>μ<sub>1</sub></th><th>μ<sub>2</sub></th><th>r<sub>A<sub>1</sub></sub></th><th>r&#773;<sub>A<sub>2</sub></sub></th><th>x<sub>1,π</sub></th><th>ε<sub>1,π</sub></th><th>s<sub>θ</sub></th><th>s<sub>A<sub>2</sub></sub></th></tr></thead><tbody><tr><td>Set 1</td><td>2.005</td><td>0.73</td><td>0.398</td><td>1.50</td><td>2.68</td><td>172.295</td><td>-0.176</td><td>0.214</td><td>0.055</td><td>0.241</td><td>32.97</td><td>5.357</td><td>15.131</td><td>4.151</td></tr><tr><td>Set 2</td><td>4.589</td><td>0.504</td><td>0.829</td><td>1.160</td><td>2.775</td><td>177.807</td><td>0.609</td><td>0.139</td><td>0.077</td><td>0.928</td><td>34.56</td><td>6.555</td><td>14.733</td><td>1.105</td></tr><tr><td>Set 3</td><td>1.795</td><td>0.535</td><td>0.398</td><td>4.537</td><td>6.097</td><td>150.576</td><td>0.176</td><td>0.678</td><td>0.032</td><td>0.29</td><td>30</td><td>4.385</td><td>17.609</td><td>6.49</td></tr></tbody></table>

## 5.1 Indicator of the invasive potential

To study the effect of denervation and the role of axons in pancreatic cancer progression, we investigate the effect of three types of denervation:

* the denervation of sympathetic axons ($A_1$),

* the denervation of sensory axons ($A_2$),