data. This can be explained by the “shift effect” from a global harmful to a protective role of the PNS, resulting from the compensation of the pro-tumor effect of denervating sympathetic axons by the anti-tumor effect of denervating sensory axons.

In addition, the model identifies different sets of parameters that may reflect biological variability of tumor innervation. These different sets highlight the importance of the level of sympathetic inhibition on tumor growth. For example, strong sympathetic inhibition on tumor growth masks the anti-tumor effect of sensory axon denervation. In this case, targeting the PNS may have little or no benefit for PDAC treatment. However, if the sympathetic inhibition of tumor growth is weak, the anti-tumor effect of the sensory denervation becomes significant. In this scenario, targeting sensory axons would be beneficial. However, it is important to note that the beneficial effects are observed with a delay. Indeed, the model highlights a latency period between denervation and observable benefits. Therefore, in terms of clinical applications for the treatment of pancreatic cancer, future knowledge of both the density of tumor innervation and the activity of sympathetic axons will be crucial for adapting the denervation strategy. This highlights the importance of a patient-specific approach to the timing of nerve block in the treatment of pancreatic cancer.

This study also highlights the fact that tumorigenesis is a constantly evolving process with inherent biological variability. To account for the biological variability and complexity of this process, the temporal components of the data need to be further investigated. On the one hand, having additional data or at least longitudinal components as data helps to reduce the uncertainty of the parameter estimation problem. On the other hand, the predictive quality of the model also depends to a large extent on the data used to calibrate it.

An improvement in the mathematical model could be to include elements of the tumor microenvironment. The complexity of pancreatic cancer tumorigenesis and its neuroregulation is also related to its interactions with the tumor microenvironment. A first step could be to incorporate spatial structure into the model. The spatio-temporal dynamics of such a model would certainly make it difficult to establish macroscopic properties and then draw biologically relevant conclusions. Conversely, a more realistic model will allow more precise investigation of biophysical properties, such as tissue stiffness and its effect on axon growth, or more complex regulatory processes, such as additional regulation by the immune system.

### CRediT authorship contribution statement

**Marie-Jose Chaaya:** Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Conceptualization. **Sophie Chauvet:** Writing – review & editing, Writing – original draft, Visualization, Methodology, Investigation, Funding acquisition, Data curation, Conceptualization. **Florence Hubert:** Writing – review & editing, Writing – original draft, Visualization, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. **Fanny Mann:** Writing – review & editing, Writing – original draft, Visualization, Investigation, Funding acquisition, Data curation, Conceptualization. **Mathieu Mezache:** Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Conceptualization. **Pierre Pudlo:** Writing – review & editing, Writing – original draft, Visualization, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

### Fundings

This research was supported by Centre National de la Recherche Scientifique (CNRS), France; Aix-Marseille University (AMU), France; grant from Institut National Du Cancer (INCa), France, Fondation Arc and Ligue contre le cancer (PAIR Pancreas, project title: “The impact of axonogenesis in pancreatic cancer”, convention number 186738) to F.M and F.H.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Appendix A. Results for the well-posedness

In this section, we give the mathematical arguments used to prove the well-posedness of the model. The proof of Theorem 1 will require four steps. First, we prove that the map $S$ is well defined, i.e the system (7) admits a unique solution. Second, we prove that the map $S$ is a contraction. 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 the following Lemmas.

**Lemma 1.** Let $T > 0$, $\theta \in C([0, T])$ such that Hypothesis 2.4 holds and $r > 0$. The Cauchy problem :

$$ egin{cases} \frac{d}{dt} y(t) = ry(t)(1 - y(t)) \left( \frac{y(t)}{\theta(t)} - 1 \right), \quad t \in [0, T], \\ y(0) = y_0 \in (0, 1), \end{cases} \qquad (22) $$

admits a unique solution $y \in C^1([0, T])$ associated to $\theta \in C([0, T])$ and $y(t) \in (0, 1), \forall t \in [0, T]$.

Moreover, assume $y_1$ (resp. $y_2$) is the unique solution of (22) associated to the Allee effect term $\theta_1 \in C([0, T])$ (resp. $\theta_2 \in C([0, T])$) where $y_1(0) = y_2(0) = y_0 \in (0, 1)$. Let $\theta_i \in C([0, T])$ such that Hypothesis 2.4 holds and $i = 1, 2$ then

$$ |y_1(t) - y_2(t)| \le C_1(\theta_2) rte^{C_2(\theta_1 - \theta_2)t} \|\theta_1 - \theta_2\|_{L^\infty([0,t])}, $$

where $C_1 > 0$ and $C_2 > 0$ only depend on $\theta_1 > 0$ and $\theta_2 > 1$.

**Proof (Lemma 1).** The local existence is a direct consequence of the Picard–Lindelöf Theorem. Moreover, since $y(t) = 0$ and $y(t) = 1$ are two stationary solutions, it implies that the solution with initial datum $y_0$ stays in$(0, 1)$ and thus is globally defined. We focus on the proof of the second item. We introduce the function $\tilde{f}: (y, \theta) \in (0, 1) \times \mathbb{R}_+^* \mapsto ry(1-y)\left(\frac{y}{\theta} - 1\right)$. Let us note that

$$ egin{aligned} |y_1(t) - y_2(t)| &\le \int_0^t \left| \tilde{f}(y_1(s), \theta_1(s)) - \tilde{f}(y_2(s), \theta_2(s)) \right| ds, \\ &\le \int_0^t \left| \tilde{f}(y_1(s), \theta_1(s)) - \tilde{f}(y_1(s), \theta_2(s)) \right| ds \\ &\quad + \int_0^t \left| \tilde{f}(y_1(s), \theta_2(s)) - \tilde{f}(y_2(s), \theta_2(s)) \right| ds. \end{aligned} $$

On the one hand, we have

$$ egin{aligned} \int_0^t \left| \tilde{f}(y_1(s), \theta_1(s)) - \tilde{f}(y_1(s), \theta_2(s)) \right| ds &\le \int_0^t r y_1(s)^2 (1 - y_1(s)) \left| \frac{1}{\theta_1(s)} - \frac{1}{\theta_2(s)} \right| ds, \\ &\le \frac{\|\theta_1 - \theta_2\|_{L^\infty([0,t])}}{\min_{\sigma \in [0,t]} |\theta_1(\sigma)\theta_2(\sigma)|} \int_0^t r y_1(s)^2 (1 - y_1(s)) ds. \end{aligned} $$

Moreover, the mapping $y \in (0, 1) \mapsto y^2(1-y)$ admits a unique maximum in $y = 2/3$. Hence, we get

$$ \int_0^t \left| \tilde{f}(y_1(s), \theta_1(s)) - \tilde{f}(y_1(s), \theta_2(s)) \right| ds \le \frac{4rt}{27 \min_{s \in [0,t]} |\theta_1(s)\theta_2(s)|} \|\theta_1 - \theta_2\|_{L^\infty([0,t])}. \qquad (23) $$

On the second hand, we have

$$ egin{aligned} \int_0^t \left| \tilde{f}(y_1(s), \theta_2(s)) - \tilde{f}(y_2(s), \theta_2(s)) \right| ds &\le \int_0^t \frac{|\tilde{f}(y_1(s), \theta_2(s)) - \tilde{f}(y_2(s), \theta_2(s))|}{|y_1(s) - y_2(s)|} |y_1(s) - y_2(s)| ds, \\ &\le \int_0^t \frac{|\tilde{f}(y_1(s), \theta_2(s)) - \tilde{f}(y_2(s), \theta_2(s))|}{|y_1(s) - y_2(s)|} |y_1(s) - y_2(s)| ds, \end{aligned} $$
• 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:

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

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 ($\mathcal{I}(T) > 0$), and conversely it is negative for an overall anti-tumoral effect of denervation ($\mathcal{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 Figure 14.

We denote by $T$ the time of observation, then the effects of the denervations are quantified by the indicator $\mathcal{I}(T)$ given by (14) 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 $\mathcal{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 $\mathcal{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 $\mathcal{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 $\mathcal{I}(T) > 0$ in the case of early denervation of $A_1$ and $A_2$.

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.

<table><caption>Table 3: Values in percentage of the invasive potential $\mathcal{I}(70)$ (cf. (14)) 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>
$$ \leq \int_{0}^{t} |\partial_{y} \bar{f}(c, \theta_{2}(s))| |y_{1}(s) - y_{2}(s)| ds, $$

where $c \in (0, 1)$ is a constant coming from the mean value theorem. Moreover, we have

$$ \begin{aligned} \left| \partial_y \bar{f}(y, \theta) \right|_{y=c} &= \left| r \left[ -1 + 2 \left( 1 + \frac{1}{\theta} \right) c - \frac{3}{\theta} c^2 \right] \right|, \\ &\leq r \max \left( \frac{(\theta - 1)^2}{3\theta} + \frac{1}{3}, 1 \right). \end{aligned} $$

Hence, we get

$$ \begin{aligned} & \int_{0}^{t} |\bar{f}(y_{1}(s), \theta_{2}(s)) - \bar{f}(y_{2}(s), \theta_{2}(s))| ds \\ &\leq \int_{0}^{t} r \max \left( \frac{(\theta_{2}(s) - 1)^{2}}{3\theta_{2}(s)} + \frac{1}{3}, 1 \right) |y_{1}(s) - y_{2}(s)| ds. \end{aligned} \quad (24) $$

Finally, using the bounds (23), (24) and Grönwall's inequality, we get the following result

$$ |y_1(t) - y_2(t)| \leq \frac{4rte^{\int_0^t r \max\left(\frac{(\theta_2(\sigma)-1)^2}{3\theta_2(\sigma)} + \frac{1}{3}, 1\right)d\sigma}}{27 \min_{s \in [0,T]} |\theta_1(s)\theta_2(s)|} \|\theta_1 - \theta_2\|_{L^\infty([0,t])}. $$

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

**Lemma 2.** Let $T > 0$, $r \in C([0, T]) \cap L^1([0, T])$ and $K \in \mathbb{R}_+^*$. The Cauchy problem :

$$ \left\{ \begin{array}{l} \frac{d}{dt} x(t) = r(t)x(t) \left(1 - \frac{x(t)}{K}\right), \quad t \in [0, T), \\ x(0) = x_0 \in (0, K), \end{array} \right. \quad (25) $$

admits a unique solution associated to $r \in C([0, T]) \cap L^1([0, T])$ given by

$$ x(t) = \frac{K}{1 + Ae^{-\int_{0}^{t} r(s) ds}}, $$

where $A = \frac{K - x_0}{x_0}$.

Moreover, assume $x_1$ (resp. $x_2$) is the unique solution of (25) associated to the growth rate $r_1 \in C([0, T]) \cap L^1([0, T])$ (resp. $r_2 \in C([0, T]) \cap L^1([0, T])$) where $r_i(t) > 0$, for $t \in [0, T]$ and $i = 1, 2$ then

$$ |x_1(t) - x_2(t)| \leq \frac{KAt}{(1 + A \min_{i=1,2} e^{-t|r_i||_{L^\infty([0,t])}})^2} \|r_1 - r_2\|_{L^\infty([0,t])}. $$

**Proof (Lemma 2).** The Picard–Lindelöf Theorem ensures that the Cauchy problem (25) admits a unique solution that stays in $(0, K)$. We focus now on the proof of the second item. In the following, we use the notation $n_i(t) = e^{-\int_0^t r_i(s) ds}$ for $i = 1, 2$ and $t \in [0, T]$. Hence, we have

$$ \begin{aligned} |x_1(t) - x_2(t)| &\leq K \frac{|A(n_1(t) - n_2(t))|}{|(1 + A n_1(t))(1 + A n_2(t))|} \\ &\leq \frac{KA}{(1 + A \min_{i=1,2} n_i(t))^2} \int_0^t |r_1(s) - r_2(s)| ds \\ &\leq \frac{KA}{(1 + A \min_{i=1,2} e^{-t|r_i||_{L^\infty([0,t])}})^2} \|r_1 - r_2\|_{L^\infty([0,t])}. \end{aligned} $$

since the following inequality holds for $(x, y) \in \mathbb{R}_+ \times \mathbb{R}_+$

$$ |e^{-x} - e^{-y}| \leq |x - y|. \quad \blacksquare $$

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

In order to be able to state the next Lemma on well-posedness of the linear PDE, we introduce the characteristics associated to (7) by :

$$ \left\{ \begin{array}{l} \frac{d}{ds} X(s, t, x; \mathcal{X}) = f(s, X(s, t, x; \mathcal{X})); \mathcal{X}), s \in \mathbb{R} \\ X(t, t, x; \mathcal{X}) = x. \end{array} \right. \quad (26) $$

**Lemma 3.** Let $\mathcal{X} \in \mathcal{P}$ (set defined in (6)) be given and assume that Hypothesis 2.1 on $Q_0$ holds. Let $f$ and $c$ defined in (7) such that the Hypotheses 2.2 and 2.3 hold. Then the PDE defined in (7) has unique solution $u \in C^1([0, T], C^1(\Omega) \cap L^1(\Omega))$ given by

$$ u(t, x; \mathcal{X}) = Q_0(X(0, t, x; \mathcal{X})) \times \exp \left( \int_0^t -c(\sigma, X(\sigma, t, x; \mathcal{X})) d\sigma \right) $$

and we have

$$ 1. \int_{\Omega} Q_0(x) dx e^{-C_1 t} \leq \|u(t)\|_{L^1(\Omega)} \leq \int_{\Omega} Q_0(x) dx e^{C_1 t}, \text{ for } t \in [0, T], $$

$$ 2. \|\partial_x u(s)\|_{L^1(\Omega)} \leq C_2 t e^{C_1 t}, \text{ for } t \in [0, T] $$

with $C_i > 0$ for $i = 1, 2$ and $C_1 = C_1^+ \sup_{0 \leq t \leq T} |h(A_1(t), A_2(t)) - N(t)|$.

Moreover, let $\mathcal{X}_i = (Q_i, A_1^i, A_2^i) \in \mathcal{P}$ and assume that $u(t, x; \mathcal{X}_i)$ is the unique solution associated to $\mathcal{X}_i$ for $i = 1, 2$, then

$$ \sup_{0 \le t \le T} \|u(t; \mathcal{X}_1) - u(t; \mathcal{X}_2)\|_{L^1(\Omega)} \le C(T)T \|\mathcal{X}_1 - \mathcal{X}_2\|_{\mathcal{P}} $$

where $0 < C(T) < \infty$.

**Proof.** First step. The existence of the unique solution of (7) under the assumptions of Lemma 3 is a classical result and can be found for instance in Chapter 6 of Perthame (2006). It relies on the existence and the regularity of the characteristics (26) that hold since $f \in C([0, T], C_1^+(\Omega))$.

Now, given $\mathcal{X} \in \mathcal{P}$, we rewrite Eq. (7) in its conservative global form (linear version of (4)):

$$ \left\{ \begin{array}{ll} \partial_t u(t, x) + \partial_x(f(t, x; \mathcal{X})u(t, x)) + g(t, x; \mathcal{X})u(t, x) = 0, & \text{for } t \in (0, T), x \in \Omega, \\ f(t, x, \mathcal{X})u(t, x) = 0, & \text{for } x \in \partial\Omega, t \in (0, T), \\ u(0, x) = Q_0(x), & \text{for } x \in \Omega. \end{array} \right. \quad (27) $$

Hence, we have that

$$ \frac{d}{dt} \int_{\Omega} u(t, x) dx = \int_{\Omega} -g(t, x; \mathcal{X})u(t, x) dx, $$

with $g(t, x; \mathcal{X})$ defined in (5). We recall that $g(t, \cdot; \mathcal{X}_i) \in C_c^1(\Omega)$, $0 \le A_i(t) \le 1$ for $t \in [0, T]$ and $i = 1, 2$ (cf. the definition of the set $\mathcal{P}$ (6)).

**Second step.** We now prove the estimates on $u$ and on $\partial_x u$. First, one can notice that $u(t, x) \ge 0$ for $t \in [0, T]$ and $x \in \Omega$ since $Q_0 \ge 0$. On the first hand, using the assumptions on $g$ from Hypothesis 2.3 and the fact that $\mathcal{X} \in \mathcal{P}$, we have

$$ \begin{aligned} \frac{d}{dt} \|u(t)\|_{L^1(\Omega)} &= \frac{d}{dt} \int_{\Omega} u(t, x) dx \leq \int_{\Omega} |g(t, x)u(t, x)| dx, \\ &\leq \|g\|_{L^{\infty}([0,T] \times \Omega)} \int_{\Omega} u(t, x) dx, \\ &\leq \|r\|_{L^{\infty}(\Omega)} \sup_{0 \le t \le T} |h(A_1(s), A_2(s)) - N(s)| \int_{\Omega} u(t, x) dx. \end{aligned} $$

On the other hand, we have for $0 \le t \le T$

$$ -\|r\|_{L^{\infty}(\Omega)} \sup_{0 \le t \le T} |h(A_1(s), A_2(s)) - N(s)| \int_{\Omega} u(t, x) dx \le \frac{d}{dt} \|u(t)\|_{L^1(\Omega)}. $$

Then, we obtain the first estimate of Lemma 3 using Gronwall's lemma. As for the second estimate, since $f(t, x) = \partial_x f(t, x) = 0$ for $x \in \partial\Omega$, we have

$$ \begin{aligned} \frac{d}{dt} \|\partial_x u(t)\|_{L^1(\Omega)} &\leq \int_{\Omega} |g(t, x)||\partial_x u(t, x)|dx + \int_{\Omega} |\partial_x g(t, x)|u(t, x)dx, \\ &\leq \|g\|_{L^{\infty}([0,T] \times \Omega)} \|\partial_x u(t)\|_{L^1(\Omega)} \\ &\quad + C(N(0), \|\partial_x g\|_{L^{\infty}([0,T] \times \Omega)}) e^{t\|g\|_{L^{\infty}([0,T] \times \Omega)}}. \end{aligned} $$

Once again, using Gronwall's lemma, we obtain the following estimate for $t \in [0, T]$

$$ \|\partial_x u(t)\|_{L^1(\Omega)} \leq C(N(0), \|\partial_x g\|_{L^{\infty}([0,T] \times \Omega)}) e^{t\|g\|_{L^{\infty}([0,T] \times \Omega)}}. $$

**Third step.** Now, we prove the stability condition. We denote $u_i := u(\cdot; \mathcal{X}_i)$ and $v := u_1 - u_2$. We obtain that $v$ satisfies the following problem

$$ \left\{ \begin{array}{ll} \partial_t v + f_1 \partial_x u_1 - f_2 \partial_x u_2 + c_1 u_1 - c_2 u_2 = 0, & \text{in } (0, T] \times \Omega, \\ v(t, x) = 0, & \text{for } x \in \partial\Omega, t \in [0, T] \\ v(0, x) = 0, & x \in \Omega, \end{array} \right. \quad (28) $$

where $f_i(t, x) = f(t, x; \mathcal{X}_i)$ and $c_i(t, x) = c(t, x; \mathcal{X}_i)$. It implies that
<figure><img src="image_23.png" /><figcaption>Figure 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>

### 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 Figure 14b, 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-posedeness 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$.
$$ \frac{d}{dt} \|v\|_{L^1(\Omega)} = \int_{\Omega} \text{sign}(u_1 - u_2) [\partial_x(f_2 u_2 - f_1 u_1) + (g_1 u_1 - g_2 u_2)] dx $$

where $g_i(t, x) = g(t, x; \chi_i)$ defined in (5).

We denote

$$ I_1(t) := \int_{\Omega} \text{sign}(u_1 - u_2) \partial_x(f_2 u_2 - f_1 u_1) dx \quad \text{and} $$

$$ I_2(t) := \int_{\Omega} \text{sign}(u_1 - u_2)(g_1 u_1 - g_2 u_2) dx. $$

On the first hand, we have

$$ \begin{aligned} I_1(t) &= \int_{\Omega} \text{sign}(u_1 - u_2) \partial_x(f_2(u_2 - u_1)) + \text{sign}(u_1 - u_2) \partial_x((f_2 - f_1)u_1) dx, \\ &\le \int_{\Omega} \partial_x(-f_2|u_2 - u_1|) dx + \int_{\Omega} |\partial_x((f_2 - f_1)u_1)| dx, \\ &\le \int_{\Omega} |\partial_x(f_2 - f_1)u_1 + (f_2 - f_1)\partial_x u_1| dx, \\ &\le \|\partial_x(f_2 - f_1)\|_{L^{\infty}(\Omega)} \|u_1\|_{L^1(\Omega)} + \|(f_2 - f_1)\|_{L^{\infty}(\Omega)} \|\partial_x u_1\|_{L^1(\Omega)}. \end{aligned} $$

Using Hypothesis 2.2 and the estimates on $u$, we obtain

$$ I_1(t) \le \|\chi_1 - \chi_2\|_p (C_0 + C_2t) e^{tC_1} \quad (29) $$

On the second hand, we have

$$ \begin{aligned} I_2(t) &\le \int_{\Omega} |g_1 - g_2| |u_1| dx + \int_{\Omega} |g_2| |u_1 - u_2| dx, \\ &\le \|u_1\|_{L^1(\Omega)} \|g_1 - g_2\|_{L^{\infty}(\Omega)} + \|g_2\|_{L^{\infty}(\Omega)} \|u_1 - u_2\|_{L^1(\Omega)}. \end{aligned} $$

Using Hypothesis 2.3 and the estimates on $u$, we obtain

$$ I_2(t) \le \|\chi_1 - \chi_2\|_p C_3 e^{tC_1} + C_1 \|v\|_{L^1(\Omega)}. \quad (30) $$

We note that the constant $C_1$ in the inequalities (29) and (30) is the same and depends on the uniform bound of the growth term $g$ which is proportional to the distance between the integral of the initial condition and the saturation constant $C(\tau_C)$.

Thanks to the inequalities (29) and (30) and using the notations of Hypotheses 2.2 and 2.3, we obtain

$$ \frac{d}{dt} \|v\|_{L^1(\Omega)} \le C_1 \|v\|_{L^1(\Omega)} + (C_0 + C_3 + C_2t) e^{tC_1} \|\chi_1 - \chi_2\|_p, $$

where

* $0 < C_0 = C(C_l(\partial_x f), \|Q_0\|_{L^1})$,

* $0 < C_1 = r_+ |C(\tau_C) - C(N(0))|,$

* $0 < C_2 = C(C_l(f), \|Q_0\|_{L^1}, \|\partial_x g\|_{L^\infty})$,

* $0 < C_3 = C(C_l(g), \|Q_0\|_{L^1})$.

Finally, using the Gronwall's Lemma on the previous inequality and taking the supremum over the interval $[0, T]$, we obtain

$$ \sup_{0 \le t \le T} \|v(s)\|_{L^1(\Omega)} \le T \left(C_0 + C_3 + \frac{C_2}{2} T\right) e^{tC_1} \|\chi_1 - \chi_2\|_p. $$

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

**Lemma 4.** Assume that $\chi := (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))$ is a solution of the system (4) and $S(\chi) = \chi$ for the mapping defined in (8). Assume Hypotheses 2.1–2.3 hold. Then

* (nonnegativity) $Q(t, x) \ge 0$ and $A_i(t) \ge 0$ for $i = 1, 2$ and $\forall (t, x) \in [0, T] \times \Omega$,

* (boundedness) $\int_{\Omega} Q_0(x) dx \le \int_{\Omega} Q(t, x) dx \le C(\tau_C)$ and $A_i(t) \le 1$ for $i = 1, 2$ and $t \in [0, T]$.

**Proof.** The nonnegativity property of the solution of system (4) is immediately obtained using the characteristics for $Q$ and noticing that $0$ is a stationary solution for $A_i$ for $i = 1, 2$. The solution can be written as

$$ Q(t, x) = Q_0(X(0, t, x)) e^{-\int_0^t c(\sigma, X(\sigma, t, x)) d\sigma} $$

and has the sign of $Q_0$. The upper bounds for $A_i$ with $i = 1, 2$ are also obtained noticing that 1 is a stationary solution of the Cauchy problems. Concerning the solution of the PDE $Q$, using Hypotheses 2.2 and 2.3, the nonnegativity of $Q$ and the fact that $\int_{\Omega} Q_0(x) \le C(\tau_C)$, we first notice that

$$ \begin{aligned} & \int_{\Omega} r(x) Q(t, x) dx(N(0) - N(t)) \le \frac{d}{dt} N(t) \\ &= \int_{\Omega} r(x) Q(t, x) dx(h(A_1(t), A_2(t)) - N(t)) \le \int_{\Omega} r(x) Q(t, x) dx(C(\tau_C) - N(t)) \end{aligned} $$

for $t \in [0, T]$ where we denote $N(t) = \int_{\Omega} Q(t, x) dx$.

Hence, we obtain the following inequalities

$$ N(0) \le N(t) \le C(\tau_C). \quad \blacksquare $$

## Appendix B. Details about the numerics

The dynamical system is implemented in Python and the algorithm can be found at (https://github.com/MarieJosec/PDE_Axons_Innerv).

### B.1. Scheme for the numerical approximation of the solution

A classical upwind scheme in space and an explicit Euler scheme in time is proposed to approximate the system that can be written as

$$ \left\{ \begin{aligned} & \partial_t Q(t, x) + \partial_x [F(x, A_1(t), A_2(t), N(t), N_c(t))Q(t, x)] \\ &= R(x, Q(t, x), A_1(t), A_2(t), N(t)), \\ & \frac{d}{dt} A_1(t) = G_1(A_1(t), N(t)), \\ & \frac{d}{dt} A_2(t) = G_2(A_2(t), N_c(t), N(t)), \end{aligned} \right. $$

Note that the function $F$ always takes non negative values.

Consider a constant step “time” discretization of the interval $[0, T]$ with $t_n = ndt$ and $n \in \{0, ..., N_T\}$, where $T = t_{N_T}$ and $dt$ is the time step. We also consider “space” discretization of the interval $[-L, L]$, $x_i = (i + \frac{1}{2})h$, and $i \in \{-p, ..., p-1\}$, $p \in \mathbb{N}^*$, where $h = \frac{L}{p}$ is the space step. We will use a finite volume approach and thus introduce the points $x_{i+\frac{1}{2}} = (i + 1)h$ and $x_{i-\frac{1}{2}} = ih$ that can be viewed as being the vertices of volume cells $M_i = [x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}]$ whose centers are the $x_i$ where $L = x_{p-\frac{1}{2}}$.

We look for an approximation $Q_i^n$ of $Q(t_n, x_i)$, $N^n$ (resp. $N_c^n$) an approximation of $\int_{-L}^L Q(t_n, x) dx$ (resp. $\int_0^L Q(t_n, x) dx$), $A_1^n$ (resp. $A_2^n$) an approximation of $A_1(t_n)$ (resp. $A_2(t_n)$). To secure the positivity of $A_1^n$ and $A_2^n$, we approximate the logarithmic function $A_1^n$ and $A_2^n$ to then go back to the approximation of $A_1^n$ and $A_2^n$ by taking the exponential function.

The scheme is then the following:

$$ \frac{h}{dt} (Q_1^{n+1} - Q_1^n) + \left( F(x_{i+\frac{1}{2}}, A_1^n, A_2^n, N^n, N_c^n) Q_i^n - F(x_{i-\frac{1}{2}}, A_1^n, A_2^n, N^n, N_c^n) Q_{i-1}^n \right) = hR(x_i, Q_i^n, A_1^n, A_2^n, N^n), $$

for $n \ge 0$ and $i \in \{-p+1, ..., p-1\}$ and

$$ \left\{ \begin{aligned} & \frac{1}{dt}(\log(A_1^{n+1}) - \log(A_1^n)) = G_1(\exp(A_1^n), N^n), \quad n \ge 0, \\ & \frac{1}{dt}(\log(A_2^{n+1}) - \log(A_2^n)) = G_2(\exp A_2^n, N_c^n), \quad n \ge 0, \\ & A_1^{n+1} = \exp(\log(A_1^{n+1})), \quad n \ge 0, \\ & A_2^{n+1} = \exp(\log(A_2^{n+1})), \quad n \ge 0. \end{aligned} \right. $$

We use a trapezoidal rule to approximate the integral terms:

$$ \left\{ \begin{aligned} N^n &= \frac{h}{2} \left[ Q_{-p}^n + \sum_{i=-p+1}^{p-2} 2Q_i^n + Q_{p-1}^n \right], \\ N_c^n &= \frac{h}{2} \left[ Q_0^n + \sum_{i=i}^{p-2} 2Q_i^n + Q_{p-1}^n \right]. \end{aligned} \right. $$

The approximation of the boundary conditions can be naturally written as
### 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 Figure 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 Figure 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 effect" is highlighted by the crossing of the green denervated and the black control curves in Figure 14.

**Outcome 2.** For both set of parameters, late denervation of both axons types inhibits tumorigenesis. In Figure 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 Figure 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 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 Figures 15 and 16.

First, we introduce the finite sequence $(s_i) \in \{5, 10, \dots, 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 Figure 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
<figure><img src="image_34.png" /><figcaption>Fig. 17. Distribution of the collection $J_3$ of parameter sets that are calibrated using our 2-dimensional criterion to fit the data and a few biological knowledge. These boxplots should be compared with the admissible ranges defined in Table 1.</figcaption></figure>

$$F(x_1, A_1^n, A_2^n, N^n, N_c^n)Q_i^n = 0, \quad n \ge 0, \quad i \in \{-p, p-1\},$$

and the initial condition is approximated by

$$\begin{cases} Q_i^0 = Q^0(x_i), & i \in \{-p, \dots, p-1\}, \\ A_1^0, & A_2^0 \text{ given.} \end{cases}$$

Note that the stability of the transport scheme requires the CFL condition

$$\frac{dt}{h} F(x_{i+\frac{1}{2}}, A_1^n, A_2^n, N^n, N_c^n) \le 1, \quad \forall n \ge, \quad \forall i = -p, \dots, p-1.$$

## B.2. Calibration through an optimization method

For a given parameter set $\theta$, the computation of our 2-dimensional criterion $G(\theta) = (G_1(\theta), G_2(\theta))$ requires two numerical simulation approximations, one acting as if we were in the control group (AA) and

a second one including a denervation treatment (OHDA). The integrals over time intervals are computed with a trapezoidal rule based on the discretization of the numerical scheme. Needless to say here that the computation of the gradient $\nabla_{\theta} G_i(\theta)$ of any coordinate of the criterion is not straightforward. This means that we have no guide to find the sets of parameters with small value of the criterion. The space we have to explore is the whole hypercube $\mathcal{H}$ of dimension 14 defined by the range column of Table 1. That is why we rely on a two-stage algorithm, see Fig. 12. At first, we explore naively the whole range of parameters, from which we adjust an instrumental Gaussian distribution truncated to $\mathcal{H}$. In a second stage, the instrumental distribution is used to draw new parameter sets that are informed by the data and the biological knowledge and we seek among those draws for the best parameter sets according to our criterion.

We start the algorithm with a massive exploration of the whole hypercube $\mathcal{H}$ of dimension 14 defined by the range column of Table 1
denervation of sympathetic axons at time $t_{A1}$ and subject to the denervation of sensory axons at time $t_{A2}$ (cf. (14) and Section 5.1). We keep track of the denervation time in the notation of the index evaluated at time $s_i$ as follows: $\mathcal{I}(s_i; t_{A1}, t_{A2})$. Hence, two kind of results are observed:

* if $\mathcal{I}(s_i; t_{A1}, t_{A2}) > 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 $\mathcal{I}(s_i; t_{A1}, t_{A2}) < 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 Figure 15 (resp. Figure 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 Figure (14d) (resp. Figure (14f)).

<figure><img src="image_24.png" /><figcaption>Figure 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>

### 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. Figure 15) and at $s_i = 40$ days for parameters set 3 (cf. Figure 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_{A1}$ 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 Figure 15 compared to Figure 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 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
<table><caption>**Table 4**
Summary statistics of the collection $J_1$ of parameter sets.</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,r</sub></th><th>ε<sub>1,r</sub></th><th>s<sub>θ</sub></th><th>s<sub>A<sub>2</sub></sub></th></tr></thead><tbody><tr><td>Mean</td><td>3.069</td><td>0.526</td><td>0.497</td><td>4.202</td><td>5.046</td><td>184.452</td><td>0.234</td><td>0.444</td><td>0.032</td><td>2.1</td><td>34.249</td><td>4.77</td><td>16.006</td><td>4.382</td></tr><tr><td>Variance</td><td>2.814</td><td>0.297</td><td>0.29</td><td>2.93</td><td>2.833</td><td>48.647</td><td>0.516</td><td>0.294</td><td>0.024</td><td>1.551</td><td>3.706</td><td>2.839</td><td>2.134</td><td>2.916</td></tr></tbody></table>

drawn as a Sobol sequence, that is to say using a quasi-Monte Carlo algorithm. This gives us the collection $J_0$ of size $n_0^{QM} = 2^{18}$ parameter sets that are distributed uniformly over the hypercube $H$. The 2-dimensional criterion is then computed for each parameter set $\theta$ in $J_0$.

Among the parameter sets in $J_0$, we select the 0.2% best sets according to the first component of the criterion $G_1(\theta)$, i.e. the 0.2% parameter sets that have the smallest value of $G_1(\theta)$. This collection is filtered again according to the second component of the criterion $G_2(\theta)$, i.e. we keep only the 20% best sets according to $G_2(\theta)$. This gives us the collection $J_1$ 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)$.

To adjust the instrumental distribution to the collection $J_1$, we compute the mean and the covariance matrix of the collection, see Table 4. The instrumental distribution is then the multivariate Gaussian distribution centered at the means given in Table 4, with a diagonal covariance matrix set according to the observed variances of this Table, and truncated (or conditioned) to stay in $H$.

The second stage of the algorithm starts with many quasi-Monte Carlo draws from the instrumental distribution, using a transformation of another Sobol sequence of dimension 14. The $n_2^{QM} = 2^{18}$ draw form the collection $J_2$ of parameter sets. The 2-dimensional criterion is then computed for each parameter set $\theta$ in $J_2$. This collection is then filtering to keep only the 0.7% best sets according to the first component of the criterion $G_1(\theta)$, and then among them the 50% best sets according to the second component of the criterion $G_2(\theta)$. This gives us the final collection $J_3$ 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)$. The final collection $J_3$ is summarized in Fig. 17.

## References

Aguirre, A.J., Bardeesy, N., Sinha, M., Lopez, L., Tuveson, D.A., Horner, J., Redston, M.S., DePinho, R.A., 2003. Activated Kras and Ink4a/Arf deficiency cooperate to produce metastatic pancreatic ductal adenocarcinoma. *Genes Dev.* 17 (24), 3112–3126.

Alexander, J., Cukierman, E., 2016. Stromal dynamic reciprocity in cancer: intricacies of fibroblastic-ECM interactions. *Curr. Opin. Cell Biol.* 42, 80–93, *Cell dynamics*.

Banh, R.S., Biancur, D.E., Yamamoto, K., Sohn, A.S., Walters, B., Kuljanin, M., Gikandi, A., Wang, H., Mancias, J.D., Schneider, R.J., et al., 2020. Neurons release serine to support mRNA translation in pancreatic cancer. *Cell* 183 (5), 1202–1218.

Biankin, A.V., Waddell, N., Kassahn, K.S., Gingras, M.-C., Muthuswamy, L.B., Johns, A.L., Miller, D.K., Wilson, P.J., Patch, A.-M., Wu, J., et al., 2012. Pancreatic cancer genomes reveal aberrations in axon guidance pathway genes. *Nature* 491 (7424), 399–405.

Chauvet, S., Hubert, F., Mann, F., Mezache, M., 2023. Tumorigenesis and axons regulation for the pancreatic cancer: a mathematical approach. *J. Theoret. Biol.* 556, 111301.

Demir, I.E., Friess, H., Ceyhan, G.O., 2015. Neural plasticity in pancreatitis and pancreatic cancer. *Nat. Rev. Gastroenterol. Hepatol.* 12 (11), 649–659.

Eftimie, R., Gibelli, L., 2020. A kinetic theory approach for modelling tumour and macrophages heterogeneity and plasticity during cancer progression. *Math. Models Methods Appl. Sci.* 30 (04), 659–683.

Guillot, J., Dominici, C., Lucchesi, A., Nguyen, H.T.T., Puget, A., Hocine, M., Rangel-Sosa, M.M., Simic, M., Nigri, J., Guillaumond, F., et al., 2022. Sympathetic axonal sprouting induces changes in macrophage populations and protects against pancreatic cancer. *Nature Commun.* 13 (1), 1985.

Klein, W.M., Hraban, R.H., Klein-Szanto, A.J., Wilentz, R.E., 2002. Direct correlation between proliferative activity and dysplasia in pancreatic intraepithelial neoplasia (panin): additional evidence for a recently proposed model of progression. *Mod. Pathol.* 15 (4), 441–447.

Liddle, R.A., 2007. The role of transient receptor potential vanilloid 1 (TRPV1) channels in pancreatitis. *Biochim. Biophys. Acta* 1772 (8), 869–878.

Lolas, G., Bianchi, A., Syrigos, K.N., 2016. Tumour-induced neoneurogenesis and perineural tumour growth: a mathematical approach. *Sci. Rep.* 6 (1), 1–10.

Perthame, B., 2006. Transport Equations in Biology. Springer Science & Business Media.

Renz, B.W., Takahashi, R., Tanaka, T., Macchini, M., Hayakawa, Y., Dantes, Z., Maurer, H.C., Chen, X., Jiang, Z., Westphalen, C.B., et al., 2018a. $\beta2$ adrenergic-neurotrophin feedforward loop promotes pancreatic cancer. *Cancer Cell* 33 (1), 75–90.

Renz, B.W., Tanaka, T., Sunagawa, M., Takahashi, R., Jiang, Z., Macchini, M., Dantes, Z., Valenti, G., White, R.A., Middelhoff, M.A., et al., 2018b. Cholinergic signaling via muscarinic receptors directly and indirectly suppresses pancreatic tumorigenesis and cancer StemnessCholinergic signaling suppresses pancreatic tumorigenesis. *Cancer Discov.* 8 (11), 1458–1473.

Saloman, J.L., Albers, K.M., Li, D., Hartman, D.J., Crawford, H.C., Muha, E.A., Rhim, A.D., Davis, B.M., 2016. Ablation of sensory neurons in a genetic model of pancreatic ductal adenocarcinoma slows initiation and progression of cancer. *Proc. Natl. Acad. Sci.* 113 (11), 3078–3083.

Sinha, S., Fu, Y.-Y., Grimont, A., Ketcham, M., Lafaro, K., Saglimbeni, J.A., Askan, G., Bailey, J.M., Melchor, J.P., Zhong, Y., et al., 2017. PanIN neuroendocrine cells promote tumorigenesis via neuronal cross-talk. *Cancer Res.* 77 (8), 1868–1879.

Winkler, F., Venkatesh, H.S., Amit, M., Batchelor, T., Demir, I.E., Deneen, B., Gutmann, D.H., Hervey-Jumper, S., Kuner, T., Mabbott, D., et al., 2023. Cancer neuroscience: state of the field, emerging directions. *Cell* 186 (8), 1689–1707.
<figure><img src="image_25.png" /><figcaption>Figure 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>

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 Figure 15 (resp. Figure 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 Figure 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 Figure 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's denervation to compensate the previous pro-tumoral effect that has taken place. In Figure 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 Figure 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 Figure 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 tumour 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 visualised 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 data. This can be explained by the "shift effect" from a global harmful to a protective role of the PNS, resulting from the compensation of the pro-tumor effect of denervating sympathetic axons by the anti-tumor effect of denervating sensory axons.

In addition, the model identifies different sets of parameters that may reflect biological variability of tumor innervation. These different sets highlight the importance of the level of sympathetic inhibition on tumor growth. For example, strong sympathetic inhibition on tumor growth masks the anti-tumor effect of sensory axon denervation. In this case, targeting the PNS may have little or no benefit for PDAC treatment. However, if the sympathetic inhibition of tumor growth is weak, the anti-tumor effect of the sensory denervation becomes significant. In this scenario, targeting sensory axons would be beneficial. However, it is important to note that the beneficial effects are observed with a delay. Indeed, the model highlights a latency period between denervation and observable benefits. Therefore, in terms of clinical applications for the treatment of pancreatic cancer, future knowledge of both the density of tumor innervation and the activity of sympathetic axons will be crucial for adapting the denervation strategy. This highlights the importance of a patient-specific approach to the timing of nerve block in the treatment of pancreatic cancer.

This study also highlights the fact that tumorigenesis is a constantly evolving process with inherent biological variability. To account for the biological variability and complexity of this process, the temporal components of the data need to be further investigated. On the one hand, having additional data or at least longitudinal components as data helps to reduce the uncertainty of the parameter estimation problem. On the other hand, the predictive quality of the model also depends to a large extent on the data used to calibrate it.

An improvement in the mathematical model could be to include elements of the tumor microenvironment. The complexity of pancreatic cancer tumorigenesis and its neuroregulation is also related to its interactions with the tumor microenvironment. A first step could be to incorporate spatial structure into the model. The spatio-temporal dynamics of such a model would certainly make it difficult to establish macroscopic properties and then draw biologically relevant conclusions. Conversely, a more realistic model will allow more precise investigation of biophysical properties, such as tissue stiffness and its effect on axon growth, or more complex regulatory processes, such as additional regulation by the immune system.

**Fundings** This research was supported by Centre National de la Recherche Scientifique (CNRS), France; Aix-Marseille University (AMU), France; grant from Institut National Du Cancer (INCa), France, Fondation Arc and Ligue contre le cancer (PAIR Pancreas, project title: "The impact of axonogenesis in pancreatic cancer", convention number 186738) to F.M and F.H.

**CRediT Conceptualization:** M.J.C, S.C, F.H, F.M, M.M, P.P; **Data curation:** S.C, F.M; **Formal analysis:** M.J.C, F.H, M.M, P.P; **Funding acquisition:** S.C, F.H, F.M, P.P; **Investigation:** M.J.C, S.C, F.H, F.M, M.M, P.P; **Methodology:** M.J.C, F.H, M.M, P.P; **Software:** M.J.C, M.M; **Visualization:** M.J.C, S.C, F.H, F.M, M.M, P.P; **Writing – original draft:** M.J.C, S.C, F.H, F.M, M.M, P.P; **Writing – review & editing:** M.J.C, S.C, F.H, F.M, M.M, P.P

# A Results for the well-posedness

In this section, we give the mathematical arguments used to prove the well-posedness of the model. The proof of Theorem 1 will require four steps. First, we prove that the map $S$ is well defined, i.e the system (7) admits a unique
solution. Second, we prove that the map $\mathcal{S}$ is a contraction. 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 the following Lemmas.

**Lemma 1.** Let $T > 0$, $\theta \in C([0, T])$ such that Hypothesis 2.4 holds and $r > 0$. The Cauchy problem :

$$ egin{cases} \frac{d}{dt}y(t) = ry(t)(1-y(t))\left(\frac{y(t)}{\theta(t)} - 1\right), & t \in [0, T], \\ y(0) = y_0 \in (0, 1), \end{cases} \qquad (15) $$

admits a unique solution $y \in C^1([0, T])$ associated to $\theta \in C([0, T])$ and $y(t) \in (0, 1)$, $\forall t \in [0, T]$.

Moreover, assume $y_1$ (resp. $y_2$) is the unique solution of (15) associated to the Allee effect term $\theta_1 \in C([0, T])$ (resp. $\theta_2 \in C([0, T])$) where $y_1(0) = y_2(0) = y_0 \in (0, 1)$. Let $\theta_i \in C([0, T])$ such that Hypothesis 2.4 holds and $i = 1, 2$ then

$$ |y_1(t) - y_2(t)| \le C_1(\theta_-)rte^{C_2(\theta_-, \theta_+)t} \|\theta_1 - \theta_2\|_{\mathcal{L}^{\infty}([0,t])}, $$

where $C_1 > 0$ and $C_2 > 0$ only depend on $\theta_-$ and $\theta_+ > 1$.

**Proof.**[Lemma 1] The local existence is a direct consequence of the Picard-Lindelöf Theorem. Moreover, since $y(t) = 0$ and $y(t) = 1$ are two stationary solutions, it implies that the solution with initial datum $y_0$ stays in $(0, 1)$ and thus is globally defined. We focus on the proof of the second item. We introduce the function $\bar{f}: (y, \theta) \in (0, 1) \times \mathbb{R}_+^* \mapsto ry(1-y)(\frac{y}{\theta}-1)$. Let us note that

$$ egin{aligned} |y_1(t) - y_2(t)| &\le \int_0^t |\bar{f}(y_1(s), \theta_1(s)) - \bar{f}(y_2(s), \theta_2(s))| ds, \\ &\le \int_0^t |\bar{f}(y_1(s), \theta_1(s)) - \bar{f}(y_1(s), \theta_2(s))| ds + \int_0^t |\bar{f}(y_1(s), \theta_2(s)) - \bar{f}(y_2(s), \theta_2(s))| ds. \end{aligned} $$

On the one hand, we have

$$ egin{aligned} \int_0^t |\bar{f}(y_1(s), \theta_1(s)) - \bar{f}(y_1(s), \theta_2(s))| ds &\le \int_0^t ry_1(s)^2(1-y_1(s)) \left| \frac{1}{\theta_1(s)} - \frac{1}{\theta_2(s)} \right| ds, \\ &\le \frac{\|\theta_1 - \theta_2\|_{\mathcal{L}^{\infty}([0,t])}}{\min_{\sigma \in [0,t]} |\theta_1(\sigma)\theta_2(\sigma)|} \int_0^t ry_1(s)^2(1-y_1(s)) ds. \end{aligned} $$

Moreover, the mapping $y \in (0, 1) \mapsto y^2(1-y)$ admits a unique maximum in $y = 2/3$. Hence, we get

$$ \int_0^t |\bar{f}(y_1(s), \theta_1(s)) - \bar{f}(y_1(s), \theta_2(s))| ds \le \frac{4rt}{27 \min_{s \in [0,t]} |\theta_1(s)\theta_2(s)|} \|\theta_1 - \theta_2\|_{\mathcal{L}^{\infty}([0,t])}. \quad (16) $$

On the second hand, we have

$$ egin{aligned} \int_0^t |\bar{f}(y_1(s), \theta_2(s)) - \bar{f}(y_2(s), \theta_2(s))| ds &\le \int_0^t \frac{|\bar{f}(y_1(s), \theta_2(s)) - \bar{f}(y_2(s), \theta_2(s))|}{|y_1(s) - y_2(s)|} |y_1(s) - y_2(s)| ds, \\ &\le \int_0^t |\partial_y \bar{f}(c, \theta_2(s))| |y_1(s) - y_2(s)| ds, \end{aligned} $$

where $c \in (0, 1)$ is a constant coming from the mean value theorem. Moreover, we have

$$ egin{aligned} |\partial_y \bar{f}(y, \theta)|_{y=c} &= |r[-1 + 2(1 + \frac{1}{\theta})c - \frac{3}{\theta}c^2]|, \\ &\le r \max \left( \frac{(\theta-1)^2}{3\theta} + \frac{1}{3}, 1 \right). \end{aligned} $$

Hence, we get

$$ \int_0^t |\bar{f}(y_1(s), \theta_2(s)) - \bar{f}(y_2(s), \theta_2(s))| ds \le \int_0^t r \max \left( \frac{(\theta_2(s) - 1)^2}{3\theta_2(s)} + \frac{1}{3}, 1 \right) |y_1(s) - y_2(s)| ds. \quad (17) $$

Finally, using the bounds (16), (17) and Grönwall's inequality, we get the following result

$$ |y_1(t) - y_2(t)| \le \frac{4rte^{\int_0^t r \max\left(\frac{(\theta_2(s)-1)^2}{3\theta_2(s)} + \frac{1}{3}, 1\right) ds}}{27 \min_{s \in [0,t]} |\theta_1(s)\theta_2(s)|} \|\theta_1 - \theta_2\|_{\mathcal{L}^{\infty}([0,t])}. $$

<figure><img src="image_26.png" /></figure>
**Lemma 2.** Let $T > 0, r \in C([0, T]) \cap \mathcal{L}^1([0, T])$ and $K \in \mathbb{R}_+^*$. The Cauchy problem :

$$ egin{cases} \frac{d}{dt}x(t) = r(t)x(t) \left(1 - \frac{x(t)}{K}\right), & t \in [0, T), \\ x(0) = x_0 \in (0, K), \end{cases} \qquad (18) $$

admits a unique solution associated to $r \in C([0, T]) \cap \mathcal{L}^1([0, T])$ given by

$$ x(t) = \frac{K}{1 + Ae^{-\int_{0}^{t} r(s)ds}}, $$

where $A = \frac{K-x_0}{x_0}$.

Moreover, assume $x_1$ (resp. $x_2$) is the unique solution of (18) associated to the growth rate $r_1 \in C([0, T]) \cap \mathcal{L}^1([0, T])$ (resp. $r_2 \in C([0, T]) \cap \mathcal{L}^1([0, T])$) where $r_i(t) > 0$, for $t \in [0, T]$ and $i = 1, 2$ then

$$ |x_1(t) - x_2(t)| \le \frac{KAt}{\left(1 + A \min_{i=1,2} e^{-t \|r_i\|_{\mathcal{L}^{\infty}([0,t])}}\right)^2} \|r_1 - r_2\|_{\mathcal{L}^{\infty}([0,t])}. $$

**Proof.** [Lemma 2] The Picard-Lindelöf Theorem ensures that the Cauchy problem (18) admits a unique solution that stays in $(0, K)$. We focus now on the proof of the second item. In the following, we use the notation $n_i(t) = e^{-\int_0^t r_i(s)ds}$ for $i = 1, 2$ and $t \in [0, T]$. Hence, we have

$$ egin{aligned} |x_1(t) - x_2(t)| &\le K \frac{|A(n_1(t) - n_2(t))|}{|(1 + An_1(t))(1 + An_2(t))|} \\ &\le \frac{KA}{\left(1 + A \min_{i=1,2} n_i(t)\right)^2} \int_0^t |r_1(s) - r_2(s)|ds \\ &\le \frac{KAt}{\left(1 + A \min_{i=1,2} e^{-t \|r_i\|_{\mathcal{L}^{\infty}([0,t])}}\right)^2} \|r_1 - r_2\|_{\mathcal{L}^{\infty}([0,t])} \end{aligned} $$

since the following inequality holds for $(x, y) \in \mathbb{R}_+ \times \mathbb{R}_+$

$$ |e^{-x} - e^{-y}| \le |x - y|. $$

In order to be able to state the next Lemma on well-posedness of the linear PDE, we introduce the characteristics associated to (7) by :

$$ \left\{ \begin{array}{l} \frac{d}{ds} X(s, t, x; \mathcal{X}) = f(s, X(s, t, x; \mathcal{X}); \mathcal{X}), \quad s \in \mathbb{R} \\ X(t, t, x; \mathcal{X}) = x. \end{array} \right. \qquad (19) $$

**Lemma 3.** Let $\mathcal{X} \in \mathcal{P}$ (set defined in (6)) be given and assume that the Hypothesis 2.1 on $Q_0$ holds. Let $f$ and $c$ defined in (7) such that the Hypotheses 2.2 and 2.3 hold. Then the PDE defined in (7) has unique solution $u \in C^1([0, T], C^1(\Omega) \cap \mathcal{L}^1(\Omega))$ given by

$$ u(t, x; \mathcal{X}) = Q_0(X(0, t, x; \mathcal{X})) \times \exp \left( \int_0^t -c(\sigma, X(\sigma, t, x; \mathcal{X})) d\sigma \right) $$

and we have

1. $\int_{\Omega} Q_0(x)dx e^{-C_1 t} \le \|u(t)\|_{\mathcal{L}^1(\Omega)} \le \int_{\Omega} Q_0(x)dx e^{C_1 t}$, for $t \in [0, T]$,

2. $\|\partial_x u(s)\|_{\mathcal{L}^1(\Omega)} \le C_2 t e^{C_1 t}$, for $t \in [0, T]$

with $C_i > 0$ for $i = 1, 2$ and $C_1 = r_+ \sup_{0 \le t \le T} |h(A_1(t), A_2(t)) - N(t)|$.

Moreover, let $\mathcal{X}_i = (Q_i, A_1^i, A_2^i) \in \mathcal{P}$ and assume that $u(t, x; \mathcal{X}_i)$ is the unique solution associated to $\mathcal{X}_i$ for $i = 1, 2$, then

$$ \sup_{0 \le t \le T} \|u(t; \mathcal{X}_1) - u(t; \mathcal{X}_2)\|_{\mathcal{L}^1(\Omega)} \le C(T)T \|\mathcal{X}_1 - \mathcal{X}_2\|_{\mathcal{P}} $$

where $0 < C(T) < \infty$.
**Proof.**

*First step.* The existence of the unique solution of (7) under the assumptions of Lemma 3 is a classical result and can be found for instance in Chapter 6 of [Perthame, 2006]. It relies on the existence and the regularity of the characteristics (19) that hold since $f \in C([0, T], C_c^1(\Omega))$.

Now, given $\mathcal{X} \in \mathcal{P}$, we rewrite the equation (7) in its conservative global form (linear version of (4)):

$$ egin{cases} 
\partial_t u(t,x) + \partial_x(f(t,x; \mathcal{X})u(t,x)) + g(t,x; \mathcal{X})u(t,x) = 0, & \text{for } t \in (0, T), x \in \Omega, \\ 
f(t,x, \mathcal{X})u(t,x) = 0, & \text{for } x \in \partial\Omega, t \in (0, T), \\ 
u(0,x) = Q_0(x), & \text{for } x \in \Omega. 
\end{cases} 
\quad (20) $$

Hence, we have that

$$ \frac{d}{dt} \int_{\Omega} u(t,x)dx = \int_{\Omega} -g(t,x; \mathcal{X})u(t,x)dx, $$

with $g(t, x; \mathcal{X})$ defined in (5). We recall that $g(t, \cdot; \mathcal{X}_i) \in C_c^1(\Omega)$, $0 \le A_i(t) \le 1$ for $t \in [0, T]$ and $i = 1, 2$ (cf. the definition of the set $\mathcal{P}$ (6)).

*Second step.* We now prove the estimates on $u$ and on $\partial_x u$. First, one can notice that $u(t, x) \ge 0$ for $t \in [0, T]$ and $x \in \Omega$ since $Q_0 \ge 0$. On the first hand, using the assumptions on $g$ from Hypothesis 2.3 and the fact that $\mathcal{X} \in \mathcal{P}$, we have

$$ \begin{aligned} \frac{d}{dt} \|u(t)\|_{L^1(\Omega)} &= \frac{d}{dt} \int_{\Omega} u(t,x)dx \le \int_{\Omega} |g(t,x)u(t,x)|dx, \\ &\le \|g\|_{L^{\infty}([0,T]\times\Omega)} \int_{\Omega} u(t,x)dx, \\ &\le \|r\|_{L^{\infty}(\Omega)} \sup_{0 \le s \le T} |h(A_1(s), A_2(s)) - N(s)| \int_{\Omega} u(t,x)dx. \end{aligned} $$

On the other hand, we have for $0 \le t \le T$

$$ -\|r\|_{L^{\infty}(\Omega)} \sup_{0 \le s \le T} |h(A_1(s), A_2(s)) - N(s)| \int_{\Omega} u(t,x)dx \le \frac{d}{dt} \|u(t)\|_{L^1(\Omega)}. $$

Then, we obtain the first estimate of Lemma 3 using Gronwall's lemma. As for the second estimate, since $f(t, x) = \partial_x f(t, x) = 0$ for $x \in \partial\Omega$, we have

$$ \begin{aligned} \frac{d}{dt} \|\partial_x u(t)\|_{L^1(\Omega)} &\le \int_{\Omega} |g(t,x)| \|\partial_x u(t,x)\| dx + \int_{\Omega} |\partial_x g(t,x)| |u(t,x)| dx, \\ &\le \|g\|_{L^{\infty}([0,T]\times\Omega)} \|\partial_x u(t)\|_{L^1(\Omega)} + C(N(0), \|\partial_x g\|_{L^{\infty}([0,T]\times\Omega)}) e^{t \|g\|_{L^{\infty}([0,T]\times\Omega)}}. \end{aligned} $$

Once again, using Gronwall's lemma, we obtain the following estimate for $t \in [0, T]$

$$ \|\partial_x u(t)\|_{L^1(\Omega)} \le C(N(0), \|\partial_x g\|_{L^{\infty}([0,T]\times\Omega)}) te^{t \|g\|_{L^{\infty}([0,T]\times\Omega)}}. $$

*Third step.* Now, we prove the stability condition. We denote $u_i := u(\cdot; \mathcal{X}_i)$ and $v := u_1 - u_2$. We obtain that $v$ satisfies the following problem

$$ egin{cases} 
\partial_t v + f_1 \partial_x u_1 - f_2 \partial_x u_2 + c_1 u_1 - c_2 u_2 = 0, & \text{in } (0, T] \times \Omega, \\ 
v(t, x) = 0, & \text{for } x \in \partial \Omega, t \in [0, T] \\ 
v(0, x) = 0, & \text{for } x \in \Omega, 
\end{cases} 
\quad (21) $$

where $f_i(t, x) = f(t, x; \mathcal{X}_i)$ and $c_i(t, x) = c(t, x; \mathcal{X}_i)$. It implies that

$$ \frac{d}{dt} \|v\|_{L^1(\Omega)} = \int_{\Omega} \text{sign}(u_1 - u_2) [\partial_x(f_2 u_2 - f_1 u_1) + (g_1 u_1 - g_2 u_2)] dx $$

where $g_i(t, x) = g(t, x; \mathcal{X}_i)$ defined in (5).

We denote

$$ I_1(t) := \int_{\Omega} \text{sign}(u_1 - u_2) \partial_x (f_2 u_2 - f_1 u_1) dx \quad \text{and} \quad I_2(t) := \int_{\Omega} \text{sign}(u_1 - u_2) (g_1 u_1 - g_2 u_2) dx. $$
On the first hand, we have

$$
\begin{aligned}
I_1(t) &= \int_{\Omega} \operatorname{sign}(u_1 - u_2) \partial_x (f_2(u_2 - u_1)) + \operatorname{sign}(u_1 - u_2) \partial_x ((f_2 - f_1)u_1) dx, \\
&\leq \int_{\Omega} \partial_x (-f_2|u_2 - u_1|) dx + \int_{\Omega} |\partial_x ((f_2 - f_1)u_1)| dx, \\
&\leq \int_{\Omega} |\partial_x (f_2 - f_1)u_1 + (f_2 - f_1)\partial_x u_1| dx, \\
&\leq \|\partial_x (f_2 - f_1)\|_{L^{\infty}(\Omega)} \|u_1\|_{L^1(\Omega)} + \|(f_2 - f_1)\|_{L^{\infty}(\Omega)} \|\partial_x u_1\|_{L^1(\Omega)}.
\end{aligned}
$$

Using Hypothesis 2.2 and the estimates on $u$, we obtain

$$I_1(t) \le \|\chi_1 - \chi_2\|_p (C_0 + C_2 t) e^{tC_1} \quad (22)$$

On the second hand, we have

$$
\begin{aligned}
I_2(t) &\le \int_{\Omega} |g_1 - g_2| |u_1| dx + \int_{\Omega} |g_2| |u_1 - u_2| dx, \\
&\le \|u_1\|_{L^1(\Omega)} \|g_1 - g_2\|_{L^{\infty}(\Omega)} + \|g_2\|_{L^{\infty}(\Omega)} \|u_1 - u_2\|_{L^1(\Omega)}.
\end{aligned}
$$

Using Hypothesis 2.3 and the estimates on $u$, we obtain

$$I_2(t) \le \|\chi_1 - \chi_2\|_p C_3 e^{tC_1} + C_1 \|v\|_{L^1(\Omega)} \quad (23)$$

We note that the constant $C_1$ in the inequalities (22) and (23) is the same and depends on the uniform bound of the growth term $g$ which is proportional to the distance between the integral of the initial condition and the saturation constant $C(\tau_C)$.

Thanks to the inequalities (22) and (23) and using the notations of Hypotheses 2.2 and 2.3, we obtain

$$\frac{d}{dt} \|v\|_{L^1(\Omega)} \le C_1 \|v\|_{L^1(\Omega)} + (C_0 + C_3 + C_2 t) e^{tC_1} \|\chi_1 - \chi_2\|_p,$$

where

* $0 < C_0 = C(C_1(\partial_x f), \|Q_0\|_{L^1})$,

* $0 < C_1 = r_+ |C(\tau_C) - C(N(0))|,$

* $0 < C_2 = C(C_1(f), \|Q_0\|_{L^1}, \|\partial_x g\|_{\infty})$,

* $0 < C_3 = C(C_1(g), \|Q_0\|_{L^1})$.

Finally, using the Gronwall's Lemma on the previous inequality and taking the supremum over the interval $[0, T]$, we obtain

$$\sup_{0 \le t \le T} \|v(s)\|_{L^1(\Omega)} \le T \left(C_0 + C_3 + \frac{C_2}{2} T\right) e^{TC_1} \|\chi_1 - \chi_2\|_p.$$

**Lemma 4.** Assume that $\chi := (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])$ is a solution of the system (4) and $S(\chi) = \chi$ for the mapping defined in (8). Assume Hypotheses 2.1,2.2, 2.3 hold. Then

* (nonnegativity) $Q(t, x) \ge 0$ and $A_i(t) \ge 0$ for $i = 1, 2$ and $\forall (t, x) \in [0, T] \times \Omega$,

* (boundedness) $\int_{\Omega} Q_0(x)dx \le \int_{\Omega} Q(t, x)dx \le C(\tau_C)$ and $A_i(t) \le 1$ for $i = 1, 2$ and $t \in [0, T]$.

**Proof.** The nonnegativity property of the solution of system (4) is immediately obtained using the characteristics for $Q$ and noticing that $0$ is a stationary solution for $A_i$ for $i = 1, 2$. The solution can be written as

$$Q(t, x) = Q_0(X(0, t, x))e^{-\int_0^t e(\sigma, X(\sigma, t, x))d\sigma}$$

and has the sign of $Q_0$. The upper bounds for $A_i$ for $i = 1, 2$ are also obtained noticing that 1 is a stationary solution of the Cauchy problems. Concerning the solution of the PDE $Q$, using Hypotheses 2.3 and 2.2, the nonnegativity of $Q$ and the fact that $\int_{\Omega} Q_0(x) \le C(\tau_C)$, we first notice that

$$\int_{\Omega} r(x)Q(t, x)dx(N(0) - N(t)) \le \frac{d}{dt}N(t) = \int_{\Omega} r(x)Q(t, x)dx(h(A_1(t), A_2(t)) - N(t)) \le \int_{\Omega} r(x)Q(t, x)dx(C(\tau_C) - N(t))$$

for $t \in [0, T]$ where we denote $N(t) = \int_{\Omega} Q(t, x)dx$.

Hence, we obtain the following inequalities

$$N(0) \le N(t) \le C(\tau_C).$$
<figure><img src="image_27.png" /></figure>

Now, we prove Theorem 1.

**Proof.**[Theorem 1] Step 1. We prove the well-posedness of the system (4). To this end, we apply a fixed point theorem for the mapping defined in (8) on the set $\mathcal{P}(T)$ defined in (6):

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

where

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

where $C(N(0))$ and $C(\tau_C)$ are positive constants and

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

The norm associated to this functional space is

$$ \|(Q, A_1, A_2)\|_{\mathcal{P}} := \left( \sup_{0 \le s \le T} \|Q(s)\|_{\mathcal{L}^1(\Omega)} + \sup_{0 \le s \le T} |A_1(s)| + \sup_{0 \le s \le T} |A_2(s)| \right). $$

For $\mathcal{X} = (Q, A_1, A_2) \in \mathcal{P}$ given, thanks to the Lemmas 1, 2 and 3, the linear problem (7) admits a unique solution

$$ (u, \tilde{A}_1, \tilde{A}_2) \in \mathcal{C}^1([0, T]; \mathcal{C}^1(\Omega) \cap \mathcal{L}^1(\Omega)) \times \mathcal{C}^1([0, T]) \times \mathcal{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}\zeta_1} \quad (24) $$

where $\zeta_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}\zeta_2} \quad (25) $$

where $\zeta_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 \mathcal{B}$. In order to prove that $\mathcal{S}$ has a fixed point, we introduce the following maps :

$$ \begin{aligned} \mathcal{F}: \mathcal{P} &\to \mathcal{B}, \quad \mathcal{F}(\mathcal{X}) = u, \\ \Gamma_1: \mathcal{P} &\to \mathcal{C}^1([0, T]), \quad \Gamma_1(\mathcal{X}) = \tilde{A}_1, \\ \Gamma_2: \mathcal{P} &\to \mathcal{C}^1([0, T]), \quad \Gamma_2(\mathcal{X}) = \tilde{A}_2. \end{aligned} $$

The contraction property of $\mathcal{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

$$ \begin{aligned} \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_-)^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_-)^2}, C_l(\theta)\right) Te^{C(\theta)T} \sup_{0 \le t \le T} \|Q_1(t) - Q_2(t)\|_{\mathcal{L}^1(\Omega)}, \end{aligned} \quad (26) $$

where $0 < C \left(r_{A_1}, \frac{1}{(\theta_-)^2}, C_l(\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

$$ \begin{aligned} \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} \left|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)\right|, \\ &\le C(A_2^0, C_l(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|, \quad (27) \\ &\le C(A_2^0, C_l(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)\|_{\mathcal{L}^1(\Omega)}, \end{aligned} $$
where $0 < C(A_2^0, C_l(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_l(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\|_{\mathcal{P}}. \qquad (28)$$

Finally, using the inequations (26), (27) and (28), we obtain that there exist two constants $0 < C_1 < \infty$ and $0 < C_2 < \infty$ such that

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

where, using the notations previously introduced and those of Hypotheses 2.1, 2.2, 2.3 and 2.4, we have

$$\bullet \ C_1(T) = C\left(C_l(r_{A_2}), A_2^0, \frac{1}{C(N(0))^2}, C(\tau_C), r_{A_1}, \frac{1}{(\theta_-)^2}, C_l(\theta), C_l(\partial_x f), \|\partial_x g\|_{\infty}, C_l(g), T\right),$$

$$\bullet \ 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_l(r_{A_2}), \frac{1}{(m_2)^2}, C(\tau_C), r_{A_1}, \frac{1}{(\theta_-)^2}, C_l(\theta), C_l(\partial_x f), \|\partial_x g\|_{\infty}, C_l(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. \qquad (29)$$

Choosing $T_1 > 0$ satisfying 24,25 and 29 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 2. 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)}. \qquad (30)$$

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

$$\mathcal{P}_1 = \mathcal{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

$$\mathcal{B}_1 := \left\{ 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] \right\}.$$

Hence, we can proceed similarly as in Step 1. Thanks to (30), 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 (24) does not cause an issue since

$$ \frac{\lambda_1}{r_+ C(\tau_c)} \le \lambda_1 \frac{\log\left(\frac{C(\tau_c)}{N(\bar{T}_1)}\right)}{\|r\|_{\infty}|C(\tau_c) - N(\bar{T}_1)|} = \frac{\lambda_1}{\|r\|_{\infty}\zeta} $$

where $\zeta \in (N(T_1); C(\tau_c))$. As for the dependance to the initial condition which appears in (25), thanks to (30), 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}\zeta} $$

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

Also, thanks to (30) and the fact that $N(0) \le N(T_1)$, the following holds

$$ 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^{-\bar{C}_2 \bar{T}}}{\bar{C}_1(\bar{T})} \le \frac{e^{-\bar{C}_2 (C(N(T_1))\bar{T})}}{C_1(\bar{T}, C(N(T_1)))}. $$

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

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

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

Step 3. 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.

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

# B Details about the numerics

The dynamical system is implemented in Python and the algorithm can be found at (https://github.com/MarieJosec/PDE_Axons_Innerv).

## B.1 Scheme for the numerical approximation of the solution

A classical upwind scheme in space and an explicit Euler scheme in time is proposed to approximate the system that can be written as

$$ \left\{ \begin{aligned} \partial_t Q(t,x) + \partial_x [F(x, A_1(t), A_2(t), N(t), N_c(t))Q(t,x)] &= R(x, Q(t,x), A_1(t), A_2(t), N(t)), \\ \frac{d}{dt} A_1(t) &= G_1(A_1(t), N(t)), \\ \frac{d}{dt} A_2(t) &= G_2(A_2(t), N_c(t), N(t)), \end{aligned} \right. $$

Note that the function $F$ always takes non negative values.

Consider a constant step "time" discretization of the interval $[0, T]$ with $t_n = ndt$ and $n \in \{0, \dots, N_T\}$, where $T = t_{N_T}$ and $dt$ is the time step. We also consider "space" discretization of the interval $[-L, L]$, $x_i = (i + \frac{1}{2})h$, and $i \in \{-p, \dots, p-1\}$, $p \in \mathbb{N}^*$, where $h = \frac{L}{p}$ is the space step. We will use a finite volume approach and thus introduce the points $x_{i+\frac{1}{2}} = (i+1)h$ and $x_{i-\frac{1}{2}} = ih$ that can be viewed as being the vertices of volume cells $M_i = [x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}]$ whose centers are the $x_i$ where $L = x_{p-\frac{1}{2}}$.

We look for an approximation $Q_i^n$ of $Q(t_n, x_i)$, $N_i^n$ (resp. $N_c^n$) an approximation of $\int_{-L}^L Q(t_n, x) dx$ (resp. $\int_0^L Q(t_n, x) dx$), $A_1^n$ (resp. $A_2^n$) an approximation of $A_1(t_n)$ (resp. $A_2(t_n)$). To secure the positivity of $A_1^n$ and $A_2^n$, we approximate the logarithmic function $A_1^n$ and $A_2^n$ to then go back to the approximation of $A_1^n$ and $A_2^n$ by taking the exponential function.

The scheme is then the following:

$$ \frac{h}{dt} (Q_i^{n+1} - Q_i^n) + \left( F(x_{i+\frac{1}{2}}, A_1^n, A_2^n, N^n, N_c^n) Q_i^n - F(x_{i-\frac{1}{2}}, A_1^n, A_2^n, N^n, N_c^n) Q_{i-1}^n \right) = hR(x_i, Q_i^n, A_1^n, A_2^n, N^n), $$
for $n \ge 0$ and $i \in \{-p+1, \dots, p-1\}$ and

$$
\left\{
\begin{array}{ll}
\frac{1}{dt}(\log(A_1^{n+1}) - \log(A_1^n)) = G_1(\exp(A_1^n), N^n), & n \ge 0, \\
\\
\frac{1}{dt}(\log(A_2^{n+1}) - \log(A_2^n)) = G_2(\exp A_2^n, N_c^n), & n \ge 0, \\
\\
A_1^{n+1} = \exp(\log(A_1^{n+1})), & n \ge 0, \\
\\
A_2^{n+1} = \exp(\log(A_2^{n+1})), & n \ge 0.
\end{array}
\right.
$$

We use a trapezoidal rule to approximate the integral terms:

$$
\left\{
\begin{array}{l}
N^n = \frac{h}{2} \left[ Q^n_{-p} + \sum_{i=-p+1}^{p-2} 2Q_i^n + Q^n_{p-1} \right], \\
\\
N_c^n = \frac{h}{2} \left[ Q_0^n + \sum_{i=i}^{p-2} 2Q_i^n + Q^n_{p-1} \right].
\end{array}
\right.
$$

The approximation of the boundary conditions can be naturally written as

$$
F(x_i, A_1^n, A_2^n, N^n, N_c^n) Q_i^n = 0, \quad n \ge 0, \quad i \in \{-p, p-1\},
$$

and the initial condition is approximated by

$$
\left\{
\begin{array}{l}
Q_i^0 = Q^0(x_i), \quad i \in \{-p, \dots, p-1\}, \\
\\
A_1^0, A_2^0 \text{ given.}
\end{array}
\right.
$$

Note that the stability of the transport scheme requires the CFL condition

$$
\frac{dt}{h} F(x_{i+\frac{1}{2}}, A_1^n, A_2^n, N^n, N_c^n) \le 1, \quad \forall n \ge, \quad \forall i = -p, \dots, p-1.
$$

B.2 Calibration through an optimization method

For a given parameter set $\theta$, the computation of our 2-dimensional criterion $G(\theta) = (G_1(\theta), G_2(\theta))$ requires two numerical simulation approximations, one acting as if we were in the control group (AA) and a second one includ-
ing a denervation treatment (OHDA). The integrals over time intervals are computed with a trapezoidal rule based on the discretization of the numerical scheme. Needless to say here that the computation of the gradient $\nabla_{\theta} G_i(\theta)$ of any coordinate of the criterion is not straightforward. This means that we have no guide to find the sets of pa-
rameters with small value of the criterion. The space we have to explore is the whole hypercube $\mathcal{H}$ of dimension
14 defined by the range column of Table 1. That is why we rely on a two-stage algorithm, see Figure 12. At first, we
explore naively the whole range of parameters, from which we adjust an instrumental Gaussian distribution trun-
cated to $\mathcal{H}$. In a second stage, the instrumental distribution is used to draw new parameter sets that are informed
by the data and the biological knowledge and we seek among those draws for the best parameter sets according to
our criterion.

We start the algorithm with a massive exploration of the whole hypercube $\mathcal{H}$ of dimension 14 defined by the range column of Table 1 drawn as a Sobol sequence, that is to say using a quasi-Monte Carlo algorithm. This gives us the collection $\mathcal{J}_0$ of size $n_0^{\text{QMC}} = 2^{18}$ parameter sets that are distributed uniformly over the hypercube $\mathcal{H}$. The 2-dimensional criterion is then computed for each parameter set $\theta$ in $\mathcal{J}_0$.

Among the parameter sets in $\mathcal{J}_0$, we select the 0.2% best sets according to the first component of the criterion $G_1(\theta)$, i.e. the 0.2% parameter sets that have the smallest value of $G_1(\theta)$. This collection is filtered again according to the second component of the criterion $G_2(\theta)$, i.e. we keep only the 20% best sets according to $G_2(\theta)$. This gives us the collection $\mathcal{J}_1$ 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)$.

To adjust the instrumental distribution to the collection $\mathcal{J}_1$, we compute the mean and the covariance matrix of the collection, see Table 4. The instrumental distribution is then the multivariate Gaussian distribution centered at the means given in Table 4, with a diagonal covariance matrix set according to the observed variances of this Table, and truncated (or conditioned) to stay in $\mathcal{H}$.

The second stage of the algorithm starts with many quasi-Monte Carlo draws from the instrumental distribu-
tion, using a transformation of another Sobol sequence of dimension 14. The $n_2^{\text{QMC}} = 2^{18}$ draw form the collection
$\mathcal{J}_2$ of parameter sets. The 2-dimensional criterion is then computed for each parameter set $\theta$ in $\mathcal{J}_2$. This collection
is then filtering to keep only the 0.7% best sets according to the first component of the criterion $G_1(\theta)$, and then
among them the 50% best sets according to the second component of the criterion $G_2(\theta)$. This gives us the final
collection $\mathcal{J}_3$ 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)$. The final collection $\mathcal{J}_3$ is summarized in Figure 17.

<table><caption>Table 4: Summary statistics of the collection $\mathcal{J}_1$ of parameter sets.</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>Mean</td><td>3.069</td><td>0.526</td><td>0.497</td><td>4.202</td><td>5.046</td><td>184.452</td><td>0.234</td><td>0.444</td><td>0.032</td><td>2.1</td><td>34.249</td><td>4.77</td><td>16.006</td><td>4.382</td></tr><tr><td>Variance</td><td>2.814</td><td>0.297</td><td>0.29</td><td>2.93</td><td>2.833</td><td>48.647</td><td>0.516</td><td>0.294</td><td>0.024</td><td>1.551</td><td>3.706</td><td>2.839</td><td>2.134</td><td>2.916</td></tr></tbody></table>

<figure><img src="image_29.png" /><figcaption>Figure 17: Distribution of the collection $\mathcal{J}_3$ of parameter sets that are calibrated using our 2-dimensional criterion to fit the data and a few biological knowledge. These boxplots should be compared with the admissible ranges defined in Table 1.</figcaption></figure>

## References

[Aguirre et al., 2003] Aguirre, A. J., Bardeesy, N., Sinha, M., Lopez, L., Tuveson, D. A., Horner, J., Redston, M. S., and DePinho, R. A. (2003). Activated kras and ink4a/arf deficiency cooperate to produce metastatic pancreatic ductal adenocarcinoma. *Genes & development*, 17(24):3112–3126.
[Alexander and Cukierman, 2016] Alexander, J. and Cukierman, E. (2016). Stromal dynamic reciprocity in cancer: intricacies of fibroblastic-ecm interactions. *Current Opinion in Cell Biology*, 42:80–93. Cell dynamics.

[Banh et al., 2020] Banh, R. S., Biancur, D. E., Yamamoto, K., Sohn, A. S., Walters, B., Kuljanin, M., Gikandi, A., Wang, H., Mancias, J. D., Schneider, R. J., et al. (2020). Neurons release serine to support mrna translation in pancreatic cancer. *Cell*, 183(5):1202–1218.

[Biankin et al., 2012] Biankin, A. V., Waddell, N., Kassahn, K. S., Gingras, M.-C., Muthuswamy, L. B., Johns, A. L., Miller, D. K., Wilson, P. J., Patch, A.-M., Wu, J., et al. (2012). Pancreatic cancer genomes reveal aberrations in axon guidance pathway genes. *Nature*, 491(7424):399–405.

[Chauvet et al., 2023] Chauvet, S., Hubert, F., Mann, F., and Mezache, M. (2023). Tumorigenesis and axons regulation for the pancreatic cancer: a mathematical approach. *Journal of Theoretical Biology*, 556:111301.

[Demir et al., 2015] Demir, I. E., Friess, H., and Ceyhan, G. O. (2015). Neural plasticity in pancreatitis and pancreatic cancer. *Nature reviews Gastroenterology & hepatology*, 12(11):649–659.

[Eftimie and Gibelli, 2020] Eftimie, R. and Gibelli, L. (2020). A kinetic theory approach for modelling tumour and macrophages heterogeneity and plasticity during cancer progression. *Mathematical Models and Methods in Applied Sciences*, 30(04):659–683.

[Guillot et al., 2022] Guillot, J., Dominici, C., Lucchesi, A., Nguyen, H. T. T., Puget, A., Hocine, M., Rangel-Sosa, M. M., Simic, M., Nigri, J., Guillaumond, F., et al. (2022). Sympathetic axonal sprouting induces changes in macrophage populations and protects against pancreatic cancer. *Nature Communications*, 13(1):1985.

[Hayakawa et al., 2017] Hayakawa, Y., Sakitani, K., Konishi, M., Asfaha, S., Niikura, R., Tomita, H., Renz, B. W., Tai-lor, Y., Macchini, M., Middelhoff, M., et al. (2017). Nerve growth factor promotes gastric tumorigenesis through aberrant cholinergic signaling. *Cancer cell*, 31(1):21–34.

[Klein et al., 2002] Klein, W. M., Hruban, R. H., Klein-Szanto, A. J., and Wilentz, R. E. (2002). Direct correlation between proliferative activity and dysplasia in pancreatic intraepithelial neoplasia (panin): additional evidence for a recently proposed model of progression. *Modern pathology*, 15(4):441–447.

[Liddle, 2007] Liddle, R. A. (2007). The role of transient receptor potential vanilloid 1 (trpv1) channels in pancreatitis. *Biochimica et Biophysica Acta (BBA)-Molecular Basis of Disease*, 1772(8):869–878.

[Lolas et al., 2016] Lolas, G., Bianchi, A., and Syrigos, K. N. (2016). Tumour-induced neoneurogenesis and perineural tumour growth: a mathematical approach. *Scientific reports*, 6(1):1–10.

[Perthame, 2006] Perthame, B. (2006). *Transport equations in biology*. Springer Science & Business Media.

[Renz et al., 2018a] Renz, B. W., Takahashi, R., Tanaka, T., Macchini, M., Hayakawa, Y., Dantes, Z., Maurer, H. C., Chen, X., Jiang, Z., Westphalen, C. B., et al. (2018a). β2 adrenergic-neurotrophin feedforward loop promotes pancreatic cancer. *Cancer cell*, 33(1):75–90.

[Renz et al., 2018b] Renz, B. W., Tanaka, T., Sunagawa, M., Takahashi, R., Jiang, Z., Macchini, M., Dantes, Z., Valenti, G., White, R. A., Middelhoff, M. A., et al. (2018b). Cholinergic signaling via muscarinic receptors directly and indirectly suppresses pancreatic tumorigenesis and cancer stemness-cholinergic signaling suppresses pancreatic tumorigenesis. *Cancer discovery*, 8(11):1458–1473.

[Saloman et al., 2016] Saloman, J. L., Albers, K. M., Li, D., Hartman, D. J., Crawford, H. C., Muha, E. A., Rhim, A. D., and Davis, B. M. (2016). Ablation of sensory neurons in a genetic model of pancreatic ductal adenocarcinoma slows initiation and progression of cancer. *Proceedings of the National Academy of Sciences*, 113(11):3078–3083.

[Sinha et al., 2017] Sinha, S., Fu, Y.-Y., Grimont, A., Ketcham, M., Lafaro, K., Saglimbeni, J. A., Askan, G., Bailey, J. M., Melchor, J. P., Zhong, Y., et al. (2017). Panin neuroendocrine cells promote tumorigenesis via neuronal cross-talk. *Cancer research*, 77(8):1868–1879.

[Winkler et al., 2023] Winkler, F., Venkatesh, H. S., Amit, M., Batchelor, T., Demir, I. E., Deneen, B., Gutmann, D. H., Hervey-Jumper, S., Kuner, T., Mabbott, D., et al. (2023). Cancer neuroscience: state of the field, emerging directions. *Cell*, 186(8):1689–1707.