# Oscillations in neuronal activity: a neuron-centered spatiotemporal model of the Unfolded Protein Response in prion diseases

Elliot M Miller, Tat Chung D. Chan, Carlos Montes-Matamoros, Omar Sharif, Laurent Pujo-Menjouet, Michael R Lindstrom

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

► To cite this version:

Elliot M Miller, Tat Chung D. Chan, Carlos Montes-Matamoros, Omar Sharif, Laurent Pujo-Menjouet, et al.. Oscillations in neuronal activity: a neuron-centered spatiotemporal model of the Unfolded Protein Response in prion diseases. Bulletin of Mathematical Biology, 2024, 86 (7), pp.82. 10.1007/s11538-024-01307-y. http://hal-04610373

HAL Id: hal-04610373

https://hal.science/hal-04610373

Submitted on 13 Jun 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
# Oscillations in neuronal activity: a neuron-centered spatiotemporal model of the Unfolded Protein Response in prion diseases

Elliot M Miller, Tat Chung D. Chan, Carlos Montes-Matamoros, Omar Sharif, Laurent Pujo-Menjouet, Michael R Lindstrom

► To cite this version:

Elliot M Miller, Tat Chung D. Chan, Carlos Montes-Matamoros, Omar Sharif, Laurent Pujo-Menjouet, et al.. Oscillations in neuronal activity: a neuron-centered spatiotemporal model of the Unfolded Protein Response in prion diseases. Bulletin of Mathematical Biology, 2024, 86 (7), pp.82. 10.1007/s11538-024-01307-y. hal-04610373

HAL Id: hal-04610373

https://hal.science/hal-04610373v1

Submitted on 13 Jun 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Oscillations in neuronal activity: a neuron-centered spatiotemporal model of the Unfolded Protein Response in prion diseases

Elliot M. Miller *

Tat Chung D. Chan †

Carlos Montes-Matamoros ‡

Omar Sharif ‡

Laurent Pujo-Menjouet §

Michael R. Lindstrom ‡

May 28, 2024

## Abstract

Many neurodegenerative diseases (NDs) are characterized by the slow spatial spread of toxic protein species in the brain. The toxic proteins can induce neuronal stress, triggering the Unfolded Protein Response (UPR), which slows or stops protein translation and can indirectly reduce the toxic load. However, the UPR may also trigger processes leading to apoptotic cell death and the UPR is implicated in the progression of several NDs. In this paper, we develop a novel mathematical model to describe the spatiotemporal dynamics of the UPR mechanism for prion diseases. Our model is centered around a single neuron, with representative proteins P (healthy) and S (toxic) interacting with heterodimer dynamics (S interacts with P to form two S's). The model takes the form of a coupled system of nonlinear reaction-diffusion equations with a delayed, nonlinear flux for P (delay from the UPR). Through the delay, we find parameter regimes that exhibit oscillations in the P- and S-protein levels. We find that oscillations are more pronounced when the S-clearance rate and S-diffusivity are small in comparison to the P-clearance rate and P-diffusivity, respectively. The oscillations become more pronounced as delays in initiating the UPR increase. We also consider quasi-realistic clinical parameters to understand how possible drug therapies can alter the course of a prion disease. We find that decreasing the production of P, decreasing the recruitment rate, increasing the diffusivity of S, increasing the UPR S-threshold, and increasing the S clearance rate appear to be the most powerful modifications to reduce the mean UPR intensity and potentially moderate the disease progression.

**Keywords:** Unfolded Protein Response, prion disease, mathematical modelling, pharmacological intervention, delay differential equations

<footnote>*College of Arts and Sciences, Culverhouse College of Business, The University of Alabama, Tuscaloosa, AL, USA</footnote>

<footnote>†Department of Mathematics, University of California Berkeley, Berkeley, CA, USA</footnote>

<footnote>‡School of Mathematical and Statistical Sciences, The University of Texas Rio Grande Valley, Edinburg, TX, USA</footnote>

<footnote>§Universite Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208, Inria, 69622 Villeurbanne, France</footnote>
# Oscillations in neuronal activity: a neuron-centered spatiotemporal model of the Unfolded Protein Response in prion diseases

Elliot M. Miller *

Tat Chung D. Chan †

Carlos Montes-Matamoros ‡

Omar Sharif ‡

Laurent Pujo-Menjouet §

Michael R. Lindstrom ‡

May 28, 2024

## Abstract

Many neurodegenerative diseases (NDs) are characterized by the slow spatial spread of toxic protein species in the brain. The toxic proteins can induce neuronal stress, triggering the Unfolded Protein Response (UPR), which slows or stops protein translation and can indirectly reduce the toxic load. However, the UPR may also trigger processes leading to apoptotic cell death and the UPR is implicated in the progression of several NDs. In this paper, we develop a novel mathematical model to describe the spatiotemporal dynamics of the UPR mechanism for prion diseases. Our model is centered around a single neuron, with representative proteins P (healthy) and S (toxic) interacting with heterodimer dynamics (S interacts with P to form two S's). The model takes the form of a coupled system of nonlinear reaction-diffusion equations with a delayed, nonlinear flux for P (delay from the UPR). Through the delay, we find parameter regimes that exhibit oscillations in the P- and S-protein levels. We find that oscillations are more pronounced when the S-clearance rate and S-diffusivity are small in comparison to the P-clearance rate and P-diffusivity, respectively. The oscillations become more pronounced as delays in initiating the UPR increase. We also consider quasi-realistic clinical parameters to understand how possible drug therapies can alter the course of a prion disease. We find that decreasing the production of P, decreasing the recruitment rate, increasing the diffusivity of S, increasing the UPR S-threshold, and increasing the S clearance rate appear to be the most powerful modifications to reduce the mean UPR intensity and potentially moderate the disease progression.

**Keywords:** Unfolded Protein Response, prion disease, mathematical modelling, pharmacological intervention, delay differential equations

<footnote>*College of Arts and Sciences, Culverhouse College of Business, The University of Alabama, Tuscaloosa, AL, USA</footnote>

<footnote>†Department of Mathematics, University of California Berkeley, Berkeley, CA, USA</footnote>

<footnote>‡School of Mathematical and Statistical Sciences, The University of Texas Rio Grande Valley, Edinburg, TX, USA</footnote>

<footnote>§Universite Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208, Inria, 69622 Villeurbanne, France</footnote>
# 1 Introduction

Neurodegenerative diseases (NDs) are devastating conditions affecting the brain and central nervous system. Two of the most common NDs are Alzheimer's Disease [1], with the gradual loss of memory and capacity to function independently, and Parkinson's Disease [2], with loss of muscle coordination and cognitive impairment. While various hypotheses exist as to the origins and etiologies of NDs, medicine is not yet able to cure most of them [3].

Despite their vastly different presentations and affected brain regions, many NDs share a common set of features. The main commonality is the presence of toxic proteins, that are thought to impair neuronal function or cause cell death, spreading through the brain [4]. During this spreading, there is a recruitment process where normal proteins can become toxic through interacting with other toxic proteins in prion diseases [5] or can become toxic through forming various oligomers [6]. In the case of AD, the toxic proteins are thought to be Amyloid-beta oligomers [6, 7]; in Parkinson's, Alpha-synuclein[8]; in prion diseases, Scrapie proteins [9]; and so on.

When subjected to the stress of a toxin, neurons are known to exhibit the Unfolded Protein Response (UPR) [10] and may temporarily shut down or limit their cellular processes, reducing their production of healthy, endogenous proteins [11]. This limiting thereby cuts down on the supply chain that could yield more toxic proteins, and allows natural clearance mechanisms of the toxins to take place. When the stress is lowered, the neuron can commence its normal cellular functions again. The activation of the UPR may also cause damage and apoptosis to the neurons. When activated, the UPR results in increasing production of the endoplasmic reticulum (ER) kinases, PERK and Ire1, which reduce protein translation. But, these proteins also lead to apoptosis through a cascade of processes [12]. Post-mortem autopsies of ND patients suggest the UPR has been activated [11]. It is also suggested that weakened UPR mechanisms may be implicated in the progression of NDs [13].

Given the innate complexities of biological systems, it is extremely difficult to study how individual mechanisms affect disease etiologies. This is where mathematical modelling can be extremely useful, giving researchers a mechanism to put together a set of assumptions and to observe the outcomes of model systems subject to those assumptions. Mathematical models have been used to study a host of different disorders and mechanisms, including AD [14, 15, 16], PD[17, 18, 19], CJD [20, 21], the UPR [22, 23, 24], and many others. The models considered cover a range of scales from studies of a single neuron or two [22], to the spatial spreading of toxic proteins in the brain through modelling the connectome and how it evolves [25].

In recent work, a simple compartment model, comprised of a nonlinear system of Delay Differential Equations (DDEs), was introduced to model the UPR in prion diseases [22]. This yielded intriguing results whereby the presence of delay was able to induce oscillations in the levels of toxic proteins, and these oscillations could be turned on/off in different parameter regimes. Moreover, the parameters could be tuned to drive the toxic protein concentration to zero.

In this paper, our objective is to extend the prior model of [22] to a delayed spatiotemporal model. We develop a nonlinear system of Reaction Diffusion Equations, with a nonlinear flux term exhibiting a delay. Our model focuses on prion diseases using heterodimer dynamics [26], centered around a single neuron. Cellular prion protein, $PrP^C$ is produced by neurons; through its interaction with misfolded, toxic scrapie prion proteins, $PrP^{Sc}$, $PrP^C$ misfolds into $PrP^{Sc}$ [27, 28]. We numerically investigate parameters that yield oscillations and which parameter modifications can reduce select measures of disease severity.

Our paper is organized as follows: we present our method of study, a mathematical model, in
# 1 Introduction

Neurodegenerative diseases (NDs) are devastating conditions affecting the brain and central nervous system. Two of the most common NDs are Alzheimer's Disease [1], with the gradual loss of memory and capacity to function independently, and Parkinson's Disease [2], with loss of muscle coordination and cognitive impairment. While various hypotheses exist as to the origins and etiologies of NDs, medicine is not yet able to cure most of them [3].

Despite their vastly different presentations and affected brain regions, many NDs share a common set of features. The main commonality is the presence of toxic proteins, that are thought to impair neuronal function or cause cell death, spreading through the brain [4]. During this spreading, there is a recruitment process where normal proteins can become toxic through interacting with other toxic proteins in prion diseases [5] or can become toxic through forming various oligomers [6]. In the case of AD, the toxic proteins are thought to be Amyloid-beta oligomers [6, 7]; in Parkinson's, Alpha-synuclein[8]; in prion diseases, Scrapie proteins [9]; and so on.

When subjected to the stress of a toxin, neurons are known to exhibit the Unfolded Protein Response (UPR) [10] and may temporarily shut down or limit their cellular processes, reducing their production of healthy, endogenous proteins [11]. This limiting thereby cuts down on the supply chain that could yield more toxic proteins, and allows natural clearance mechanisms of the toxins to take place. When the stress is lowered, the neuron can commence its normal cellular functions again. The activation of the UPR may also cause damage and apoptosis to the neurons. When activated, the UPR results in increasing production of the endoplasmic reticulum (ER) kinases, PERK and Ire1, which reduce protein translation. But, these proteins also lead to apoptosis through a cascade of processes [12]. Post-mortem autopsies of ND patients suggest the UPR has been activated [11]. It is also suggested that weakened UPR mechanisms may be implicated in the progression of NDs [13].

Given the innate complexities of biological systems, it is extremely difficult to study how individual mechanisms affect disease etiologies. This is where mathematical modelling can be extremely useful, giving researchers a mechanism to put together a set of assumptions and to observe the outcomes of model systems subject to those assumptions. Mathematical models have been used to study a host of different disorders and mechanisms, including AD [14, 15, 16], PD[17, 18, 19], CJD [20, 21], the UPR [22, 23, 24], and many others. The models considered cover a range of scales from studies of a single neuron or two [22], to the spatial spreading of toxic proteins in the brain through modelling the connectome and how it evolves [25].

In recent work, a simple compartment model, comprised of a nonlinear system of Delay Differential Equations (DDEs), was introduced to model the UPR in prion diseases [22]. This yielded intriguing results whereby the presence of delay was able to induce oscillations in the levels of toxic proteins, and these oscillations could be turned on/off in different parameter regimes. Moreover, the parameters could be tuned to drive the toxic protein concentration to zero.

In this paper, our objective is to extend the prior model of [22] to a delayed spatiotemporal model. We develop a nonlinear system of Reaction Diffusion Equations, with a nonlinear flux term exhibiting a delay. Our model focuses on prion diseases using heterodimer dynamics [26], centered around a single neuron. Cellular prion protein, $PrP^C$ is produced by neurons; through its interaction with misfolded, toxic scrapie prion proteins, $PrP^{Sc}$, $PrP^C$ misfolds into $PrP^{Sc}$ [27, 28]. We numerically investigate parameters that yield oscillations and which parameter modifications can reduce select measures of disease severity.

Our paper is organized as follows: we present our method of study, a mathematical model, in
Section 2; an exploration of the model and its general behavior is covered in Section 3; we focus on the biological parameters and treatment implications for prion diseases in Section 4; and our work is concluded in Section 5. The Appendices A and B provide additional mathematical details not located in the main text.

# 2 Materials and Methods

We focus on gaining an understanding into the UPR mechanism and its relation to prion diseases at the neuron scale. In isolation, different mechanisms have been studied experimentally; however, at present, an in-depth knowledge of *in vivo* parameters of biological significance is lacking. For this reason, we combine multiple well-established biological processes into a *mathematical model*, to simulate a system and glean understanding into its dynamics. The hope is that enough of the important biology is present in the model for it to provide clinically relevant insights, even if only in approximation. Our experiments then take the form of numerical simulations, where each parameter can be carefully controlled.

A reader less focused on the mathematics, with more interest in the biology, should consider reading Section 2.1.1 and then reading Sections 3 and 4. It would also be useful to refer to Table 1 to see the different parameters studied in our model.

## 2.1 Model

This section concerns developing a mathematical model.

### 2.1.1 Overview

Our objectives with this model are to uncover the spatiotemporal dynamics of the Unfolded Protein Response. For simplicity, we consider two representative protein species, P (healthy proteins, different forms of the cellular prion protein, $PrP^C$) and S (misfolded proteins, different forms of the toxic scrapie protein, $PrP^{Sc}$). The S-proteins can recruit P-proteins to become S-proteins. We do not model higher-order structures — dimers, higher-order oligomers, nuclei, and fibrils are not present. We argue that the capacity for S to recruit P is at least representative of the more detailed biochemistry whereby misfolded proteins can form oligomers and nuclei, which can fragment to generate further recruitment. This is often referred to as a heterodimer model.

Our model centers around one neuron, with dynamics taking place in the interstitial fluid surrounding it. The P-proteins are produced within the neuron and released through the membrane into the intercellular space. We note cellular prion protein and scrapie protein can be found as both membrane-bound and extracellular forms [29, 30], but we focus upon the latter here to build a simple model. Through a buildup of S near the membrane, the P-production is lowered (due to the UPR) after some delay representing the processing time needed to stop/slow translation.

Figure 1 depicts the relevant mechanisms. Our simulations center around the geometries depicted in Figure 2.

### 2.1.2 Derivation

We denote *P* the concentration of P-proteins and *S* the concentration of S-proteins. We use *x* for position and *t* for time.
Section 2; an exploration of the model and its general behavior is covered in Section 3; we focus on the biological parameters and treatment implications for prion diseases in Section 4; and our work is concluded in Section 5. The Appendices A and B provide additional mathematical details not located in the main text.

# 2 Materials and Methods

We focus on gaining an understanding into the UPR mechanism and its relation to prion diseases at the neuron scale. In isolation, different mechanisms have been studied experimentally; however, at present, an in-depth knowledge of *in vivo* parameters of biological significance is lacking. For this reason, we combine multiple well-established biological processes into a *mathematical model*, to simulate a system and glean understanding into its dynamics. The hope is that enough of the important biology is present in the model for it to provide clinically relevant insights, even if only in approximation. Our experiments then take the form of numerical simulations, where each parameter can be carefully controlled.

A reader less focused on the mathematics, with more interest in the biology, should consider reading Section 2.1.1 and then reading Sections 3 and 4. It would also be useful to refer to Table 1 to see the different parameters studied in our model.

## 2.1 Model

This section concerns developing a mathematical model.

### 2.1.1 Overview

Our objectives with this model are to uncover the spatiotemporal dynamics of the Unfolded Protein Response. For simplicity, we consider two representative protein species, P (healthy proteins, different forms of the cellular prion protein, $PrP^C$) and S (misfolded proteins, different forms of the toxic scrapie protein, $PrP^{Sc}$). The S-proteins can recruit P-proteins to become S-proteins. We do not model higher-order structures — dimers, higher-order oligomers, nuclei, and fibrils are not present. We argue that the capacity for S to recruit P is at least representative of the more detailed biochemistry whereby misfolded proteins can form oligomers and nuclei, which can fragment to generate further recruitment. This is often referred to as a heterodimer model.

Our model centers around one neuron, with dynamics taking place in the interstitial fluid surrounding it. The P-proteins are produced within the neuron and released through the membrane into the intercellular space. We note cellular prion protein and scrapie protein can be found as both membrane-bound and extracellular forms [29, 30], but we focus upon the latter here to build a simple model. Through a buildup of S near the membrane, the P-production is lowered (due to the UPR) after some delay representing the processing time needed to stop/slow translation.

Figure 1 depicts the relevant mechanisms. Our simulations center around the geometries depicted in Figure 2.

### 2.1.2 Derivation

We denote *P* the concentration of P-proteins and *S* the concentration of S-proteins. We use *x* for position and *t* for time.
From a chemical reaction perspective, we assume that

• P and S combine at a rate $c \ge 0$ to form two separate S proteins;

• P is cleared at a rate $a \ge 0$;

• S is cleared at a rate $b \ge 0$, with $b \le a$;

• spontaneously, P can misfold to become S at a rate $f \ge 0$; and

• P and S have diffusivities within the interneuronal space of $D_P$ and $D_S$, respectively, with $0 \le D_S \le D_P$.

The fact that $D_S \le D_P$ stems from the fact that we expect, in general, a single S-protein is at least as massive as a single P-protein (recall we are not directly modelling oligomers but S represents misfolded proteins of all sizes). We also anticipate that $b \le a$ as larger proteins likely take longer to break down. Given the relative rarity of prion diseases, we expect that $f \ll a$.

To account for the UPR mechanism, we prescribe that the rate that P-proteins are produced and released decreases as the S concentration at the cell membrane increases. More precisely, the magnitude of the flux, $J$, of P from the neuron on the neuron boundary at time $t$ is given by the proportionality relation

$$J \propto \left(1 + \left(\frac{\langle S(\cdot, t - t_d) \rangle}{S_c}\right)^m\right)^{-1}, \quad (1)$$

with the constant of proportionality $A > 0$ being the maximum possible P-flux and where $S_c > 0$, $m > 0$, and $t_d \ge 0$ are prescribed constants. The $\langle \cdot \rangle$ denotes the mean value over the membrane. The value $S_c$ is more or less a sensitivity of the neurons to S: when $\langle S \rangle > S_c$, the flux may decrease very rapidly and when $\langle S \rangle < S_c$, the flux stays near its maximum. The value $m$ controls how rapidly the neuron switches from maximum P-flux to near zero P-flux. The delay $t_d$ models the fact that during stress, the P-production cannot shut down immediately; a cascade of signals needs to be transmitted to decrease P-production, which is modelled by $t_d$. The same principle applies to increasing P production when S is cleared.

We consider the system in $d$-dimensional space. While $d = 3$ is most natural, certain geometric arrangements of neurons could be modelled as being 1 or 2 dimensional. We let the neuron occupy a bounded, closed, and connected region $\Theta \subset \mathbb{R}^d$ containing the origin. In $d = 1$ we assume, without loss of generality, that $\Theta = [-R, R]$ for some length $R > 0$. For $d = 2, 3$, we further assume a smooth boundary and that a characteristic length scale for $\Theta$ is $R > 0$. The computational domain is the region $\Omega := \mathbb{R}^d \setminus \Theta$ and the cell membrane is given by $\partial\Theta = \partial\Omega$. With $d = 1$, we only model the region $(R, \infty)$ so that $\langle S(\cdot, t - t_d) \rangle = S(R, t - t_d)$. For $d = 2, 3$, $\langle S(\cdot, t - t_d) \rangle = \frac{1}{|\partial\Omega|} \int_{\partial\Omega} S(x, t - t_d)dx$.

For boundary conditions, we assume that P and S tend to 0 as $|x| \to \infty$ (exponential decay is common in diffusion problems) and that there is no flux of S on the boundary $\partial\Omega$. The flux of P on $\partial\Omega$ is based on the UPR assumptions outlined above.

From our preceding assumptions, our model amounts to the following system on $\Omega \times (0, \infty)$:

$$P_{,t} = \overbrace{D_P \triangle P}^{\text{diffusion}} - \overbrace{cPS}^{\text{recruitment loss}} - \overbrace{fP}^{\text{spontaneous misfolding}} - \overbrace{aP}^{\text{breakdown}} \quad (2)$$

$$S_{,t} = \overbrace{D_S \triangle S}^{\text{diffusion}} + \overbrace{cPS}^{\text{recruitment gain}} + \overbrace{fP}^{\text{spontaneous misfolding}} - \overbrace{bS}^{\text{breakdown}} \quad (3)$$
From a chemical reaction perspective, we assume that

• P and S combine at a rate $c \ge 0$ to form two separate S proteins;

• P is cleared at a rate $a \ge 0$;

• S is cleared at a rate $b \ge 0$, with $b \le a$;

• spontaneously, P can misfold to become S at a rate $f \ge 0$; and

• P and S have diffusivities within the interneuronal space of $D_P$ and $D_S$, respectively, with $0 \le D_S \le D_P$.

The fact that $D_S \le D_P$ stems from the fact that we expect, in general, a single S-protein is at least as massive as a single P-protein (recall we are not directly modelling oligomers but S represents misfolded proteins of all sizes). We also anticipate that $b \le a$ as larger proteins likely take longer to break down. Given the relative rarity of prion diseases, we expect that $f \ll a$.

To account for the UPR mechanism, we prescribe that the rate that P-proteins are produced and released decreases as the S concentration at the cell membrane increases. More precisely, the magnitude of the flux, $J$, of P from the neuron on the neuron boundary at time $t$ is given by the proportionality relation

$$J \propto \left(1 + \left(\frac{\langle S(\cdot, t-t_d) \rangle}{S_c}\right)^m\right)^{-1}, \quad (1)$$

with the constant of proportionality $A > 0$ being the maximum possible P-flux and where $S_c > 0$, $m > 0$, and $t_d \ge 0$ are prescribed constants. The $\langle \cdot \rangle$ denotes the mean value over the membrane. The value $S_c$ is more or less a sensitivity of the neurons to S: when $\langle S \rangle > S_c$, the flux may decrease very rapidly and when $\langle S \rangle < S_c$, the flux stays near its maximum. The value $m$ controls how rapidly the neuron switches from maximum P-flux to near zero P-flux. The delay $t_d$ models the fact that during stress, the P-production cannot shut down immediately; a cascade of signals needs to be transmitted to decrease P-production, which is modelled by $t_d$. The same principle applies to increasing P production when S is cleared.

We consider the system in $d$-dimensional space. While $d = 3$ is most natural, certain geometric arrangements of neurons could be modelled as being 1 or 2 dimensional. We let the neuron occupy a bounded, closed, and connected region $\Theta \subset \mathbb{R}^d$ containing the origin. In $d = 1$ we assume, without loss of generality, that $\Theta = [-R, R]$ for some length $R > 0$. For $d = 2, 3$, we further assume a smooth boundary and that a characteristic length scale for $\Theta$ is $R > 0$. The computational domain is the region $\Omega := \mathbb{R}^d \setminus \Theta$ and the cell membrane is given by $\partial\Theta = \partial\Omega$. With $d = 1$, we only model the region $(R, \infty)$ so that $\langle S(\cdot, t-t_d) \rangle = S(R, t-t_d)$. For $d = 2, 3$, $\langle S(\cdot, t-t_d) \rangle = \frac{1}{|\partial\Omega|} \int_{\partial\Omega} S(x, t-t_d)dx$.

For boundary conditions, we assume that P and S tend to 0 as $|x| \to \infty$ (exponential decay is common in diffusion problems) and that there is no flux of S on the boundary $\partial\Omega$. The flux of P on $\partial\Omega$ is based on the UPR assumptions outlined above.

From our preceding assumptions, our model amounts to the following system on $\Omega \times (0, \infty)$:

$$P_{,t} = \overbrace{D_P \triangle P}^{\text{diffusion}} - \overbrace{cPS}^{\text{recruitment loss}} - \overbrace{fP}^{\text{spontaneous misfolding}} - \overbrace{aP}^{\text{breakdown}} \quad (2)$$

$$S_{,t} = \overbrace{D_S \triangle S}^{\text{diffusion}} + \overbrace{cPS}^{\text{recruitment gain}} + \overbrace{fP}^{\text{spontaneous misfolding}} - \overbrace{bS}^{\text{breakdown}} \quad (3)$$
<figure><img src="image_3.png" /><figcaption>Figure 1: Representative sketch of model. Left: a cell produces protein P that is released into the intercellular space and diffuses with a rate $D_P$. P may be recruited to become S that diffuses with a rate $D_S$. The neuron is sensitive to S at a concentration $S_c$ and will reduce its maximum flux of P from A in the presence of S. Right: Spontaneous misfolding from P- to S-proteins occurs at a rate f, recruitment by interaction with an S-protein happens with a rate c, and degradation rates for P- and S-proteins are respectively a and b.</figcaption></figure>

with

$$P, S \to 0, \quad |x| \to \infty, \qquad (4)$$

$$-D_P \nabla P(x, t) \cdot \hat{n} = \frac{-A}{1 + \left( \frac{S(x, t - t_d)}{S_c} \right)^m}, \quad x \in \partial\Omega, t \ge 0 \qquad (5)$$

$$-D_S \nabla S \cdot \hat{n} = 0, \quad \text{on } \partial\Omega, \qquad (6)$$

$$(P(x, 0), S(x, 0)) = (P_0(x), S_0(x)), \quad x \in \bar{\Omega}, \qquad (7)$$

$$S(x, t) = S_{-}(x, t), \quad x \in \partial\Omega, t \in [-t_d, 0]. \qquad (8)$$

Note that $\hat{n}$, the outward unit normal is based on the frame of $\Omega$. Hence, for the influx, we have $-A$ in Eq. 5. We assume that $S_0(x) = S_-(x, 0)$ on $\partial\Omega$. We prescribe that $P_0$ and $S_0$ satisfy the proper boundary conditions, namely $P_0, S_0 \to 0$ as $|x| \to \infty$, $\nabla S_0 \cdot \hat{n} = 0$ on $\partial\Omega$, and $-D_P \nabla P_0(x) \cdot \hat{n} = \frac{-A}{1 + \left( \frac{S_-(x, t - t_d)}{S_c} \right)^m}$ on $\partial\Omega$. We note this is likely not strictly necessary as equations with diffusion tend to smooth out irregularities instantaneously. We also assume a finite total amount of protein in the system so that $\int_{\Omega} P_0(x)dx, \int_{\Omega} S_0(x)dx < \infty$. See Table 1 for a listing of all dimensional parameters. Figure 1 depicts this system.

Since a sustained UPR mechanism itself may cause damage to the cells [31, 32], we are interested in coming up with a proxy for what fraction of time the mechanism is active and to what extent. To that end, given the P-flux magnitude, $J$, we note that the quantity

$$h = 1 - J/A \in [0, 1] \qquad (9)$$

effectively describes whether the UPR is on or off at any given time and the closer $h$ is to 1, the more strongly the UPR is activated. Then, the time average $h_{\text{avg}}$, gives a loose estimate for the fraction of time the UPR is active — an effective UPR intensity.
<figure><img src="image_2.png" /><figcaption>Figure 1: Representative sketch of model. Left: a cell produces protein P that is released into the intercellular space and diffuses with a rate $D_P$. P may be recruited to become S that diffuses with a rate $D_S$. The neuron is sensitive to S at a concentration $S_c$ and will reduce its maximum flux of P from A in the presence of S. Right: Spontaneous misfolding from P- to S-proteins occurs at a rate f, recruitment by interaction with an S-protein happens with a rate c, and degradation rates for P- and S-proteins are respectively a and b.</figcaption></figure>

with

$$P, S \to 0, \quad |x| \to \infty, \qquad (4)$$

$$-D_P \nabla P(x, t) \cdot \hat{n} = \frac{-A}{1 + \left( \frac{S(\cdot, t-t_d)}{S_c} \right)^m}, \quad x \in \partial\Omega, t \ge 0 \qquad (5)$$

$$-D_S \nabla S \cdot \hat{n} = 0, \quad \text{on } \partial\Omega, \qquad (6)$$

$$(P(x, 0), S(x, 0)) = (P_0(x), S_0(x)), \quad x \in \bar{\Omega}, \qquad (7)$$

$$S(x, t) = S_{-}(x, t), \quad x \in \partial\Omega, t \in [-t_d, 0]. \qquad (8)$$

Note that $\hat{n}$, the outward unit normal is based on the frame of $\Omega$. Hence, for the influx, we have $-A$ in Eq. 5. We assume that $S_0(x) = S_-(x, 0)$ on $\partial\Omega$. We prescribe that $P_0$ and $S_0$ satisfy the proper boundary conditions, namely $P_0, S_0 \to 0$ as $|x| \to \infty$, $\nabla S_0 \cdot \hat{n} = 0$ on $\partial\Omega$, and $-D_P \nabla P_0(x) \cdot \hat{n} = \frac{-A}{1+\left(\frac{S_-(x, t-t_d)}{S_c}\right)^m}$ on $\partial\Omega$. We note this is likely not strictly necessary as equations with diffusion tend to smooth out irregularities instantaneously. We also assume a finite total amount of protein in the system so that $\int_{\Omega} P_0(x)dx, \int_{\Omega} S_0(x)dx < \infty$. See Table 1 for a listing of all dimensional parameters. Figure 1 depicts this system.

Since a sustained UPR mechanism itself may cause damage to the cells [31, 32], we are interested in coming up with a proxy for what fraction of time the mechanism is active and to what extent. To that end, given the P-flux magnitude, $J$, we note that the quantity

$$h = 1 - J/A \in [0, 1] \qquad (9)$$

effectively describes whether the UPR is on or off at any given time and the closer $h$ is to 1, the more strongly the UPR is activated. Then, the time average $h_{\text{avg}}$, gives a loose estimate for the fraction of time the UPR is active — an effective UPR intensity.
<table><caption>Table 1: All dimensional parameters in the model where M = mass, L = length, and T = time.</caption><thead><tr><th>Parameter</th><th>SI Base Units</th><th>Interpretation</th></tr></thead><tbody><tr><td>D<sub>P</sub></td><td>L<sup>2</sup>/T</td><td>diffusion coefficient of P</td></tr><tr><td>D<sub>S</sub></td><td>L<sup>2</sup>/T</td><td>diffusion coefficient of S</td></tr><tr><td>c</td><td>L<sup>3</sup>/(MT)</td><td>rate of combination of P and S</td></tr><tr><td>f</td><td>1/T</td><td>spontaneous misfolding rate of P</td></tr><tr><td>a</td><td>1/T</td><td>clearance rate of P</td></tr><tr><td>b</td><td>1/T</td><td>clearance rate of S</td></tr><tr><td>t<sub>d</sub></td><td>T</td><td>delay time of the UPR</td></tr><tr><td>A</td><td>M/(L<sup>2</sup>T)</td><td>maximum possible P-flux</td></tr><tr><td>S<sub>c</sub></td><td>M/L<sup>3</sup></td><td>critical S-concentration</td></tr><tr><td>R</td><td>L</td><td>representative neuron size/radius</td></tr></tbody></table>

<table><caption>Table 2: Scales chosen where M = mass, L = length, and T = time.</caption><thead><tr><th>Scale</th><th>SI Base Units</th><th>Definition</th><th>Interpretation</th></tr></thead><tbody><tr><td>t̄</td><td>T</td><td>1/a</td><td>characteristic time for P to be cleared</td></tr><tr><td>x̄</td><td>L</td><td>√D<sub>p</sub>/a</td><td>characteristic distance P diffuses before being cleared</td></tr><tr><td>P</td><td>M/L<sup>3</sup></td><td>A/√D<sub>p</sub>a</td><td>characteristic scale of P as measured on the diffusive scale x̄, given the flux</td></tr><tr><td>S</td><td>M/L<sup>3</sup></td><td>S<sub>c</sub></td><td>critical S-concentration to trigger decreased P-production</td></tr></tbody></table>

An alternative point of study would be to consider the frequency of UPR activations over times on (0, T). We define

$$W_T = \frac{1}{T} \left| \{t|t \text{ is a strict local maximum for } h \text{ on } (0,T)\} \right|, \quad (10)$$

which computes an average frequency of h-peaks on (0, T).

### 2.1.3 Nondimensionalization

Nondimensionalizing and rescaling our system via $x = \bar{x}\tilde{x}$, $\Omega = \bar{x}\tilde{\Omega}$, $t = \bar{t}\tilde{t}$, $P(x,t) = \bar{P}p(\tilde{x}, \tilde{t})$, and $S(x,t) = \bar{S}s(\tilde{x}, \tilde{t})$ (see Table 2), after simplifying and removing tildes, on $\Omega \times (0, \infty)$, we have

$$p_{,t} = \Delta p - \gamma p s - (1 + \sigma)p \quad (11)$$

$$s_{,t} = \delta \Delta s + \eta \gamma p s + \eta \sigma p - \beta s \quad (12)$$
<table><caption>Table 1: All dimensional parameters in the model where M = mass, L = length, and T = time.</caption><thead><tr><th>Parameter</th><th>SI Base Units</th><th>Interpretation</th></tr></thead><tbody><tr><td>D<sub>P</sub></td><td>L<sup>2</sup>/T</td><td>diffusion coefficient of P</td></tr><tr><td>D<sub>S</sub></td><td>L<sup>2</sup>/T</td><td>diffusion coefficient of S</td></tr><tr><td>c</td><td>L<sup>3</sup>/(MT)</td><td>rate of combination of P and S</td></tr><tr><td>f</td><td>1/T</td><td>spontaneous misfolding rate of P</td></tr><tr><td>a</td><td>1/T</td><td>clearance rate of P</td></tr><tr><td>b</td><td>1/T</td><td>clearance rate of S</td></tr><tr><td>t<sub>d</sub></td><td>T</td><td>delay time of the UPR</td></tr><tr><td>A</td><td>M/(L<sup>2</sup>T)</td><td>maximum possible P-flux</td></tr><tr><td>S<sub>c</sub></td><td>M/L<sup>3</sup></td><td>critical S-concentration</td></tr><tr><td>R</td><td>L</td><td>representative neuron size/radius</td></tr></tbody></table>

<table><caption>Table 2: Scales chosen where M = mass, L = length, and T = time.</caption><thead><tr><th>Scale</th><th>SI Base Units</th><th>Definition</th><th>Interpretation</th></tr></thead><tbody><tr><td>t̄</td><td>T</td><td>1/a</td><td>characteristic time for P to be cleared</td></tr><tr><td>x̄</td><td>L</td><td>√D<sub>p</sub>/a</td><td>characteristic distance P diffuses before being cleared</td></tr><tr><td>P</td><td>M/L<sup>3</sup></td><td>A/√D<sub>p</sub>a</td><td>characteristic scale of P as measured on the diffusive scale x̄, given the flux</td></tr><tr><td>S</td><td>M/L<sup>3</sup></td><td>S<sub>c</sub></td><td>critical S-concentration to trigger decreased P-production</td></tr></tbody></table>

An alternative point of study would be to consider the frequency of UPR activations over times on (0, T). We define

$$W_T = \frac{1}{T} |\{t| t \text{ is a strict local maximum for } h \text{ on } (0, T)\}|, \quad (10)$$

which computes an average frequency of h-peaks on (0, T).

### 2.1.3 Nondimensionalization

Nondimensionalizing and rescaling our system via $x = \bar{x}\tilde{x}$, $\Omega = \bar{x}\tilde{\Omega}$, $t = \bar{t}\tilde{t}$, $P(x, t) = \bar{P}p(\tilde{x}, \tilde{t})$, and $S(x, t) = \bar{S}s(\tilde{x}, \tilde{t})$ (see Table 2), after simplifying and removing tildes, on $\Omega \times (0, \infty)$, we have

$$p_{,t} = \Delta p - \gamma p s - (1 + \sigma)p \quad (11)$$

$$s_{,t} = \delta \Delta s + \eta \gamma p s + \eta \sigma p - \beta s \quad (12)$$
<table><caption>Table 3: Dimensionless parameters appearing in the model.</caption><thead><tr><th>Parameter</th><th>Definition</th><th>Interpretation</th></tr></thead><tbody><tr><td>γ</td><td>cS<sub>c</sub>/a</td><td>ratio of (the rate P is converted to S when S concentration is S<sub>c</sub>) to (the rate P is cleared)</td></tr><tr><td>σ &ll; 1</td><td>f/a</td><td>ratio of the rate P misfolds to the rate it is cleared</td></tr><tr><td>δ &le; 1</td><td>D<sub>S</sub>/D<sub>P</sub></td><td>ratio of S diffusivity to P diffusivity</td></tr><tr><td>η</td><td>A/S<sub>c</sub>√D<sub>P</sub>a</td><td>ratio of characteristic P-concentration to S<sub>c</sub></td></tr><tr><td>β &le; 1</td><td>b/a</td><td>ratio of rate S is cleared to the rate P is cleared</td></tr><tr><td>τ</td><td>at<sub>d</sub></td><td>delay relative to P-clearance time</td></tr><tr><td>ρ</td><td>R√a/√D<sub>P</sub></td><td>ratio of (the characteristic cell size) to (the characteristic distance P travels before being cleared)</td></tr><tr><td>m &gt; 0</td><td>chosen value</td><td>controls abruptness of switch</td></tr></tbody></table>

with

$$p, s \to 0, \quad |x| \to \infty \qquad (13)$$

$$-\nabla p \cdot \hat{n} = \frac{-1}{1 + (s(\cdot, t - \tau))^m} \quad x \in \partial\Omega, t \ge 0 \qquad (14)$$

$$-\nabla s \cdot \hat{n} = 0, \quad \text{on } \partial\Omega \qquad (15)$$

$$(p(x, 0), s(x, 0)) = (p_{0}^{*}(x), s_{0}^{*}(x)), \quad x \in \bar{\Omega} \qquad (16)$$

$$s(x, t) = s_{-}^{*}(x, t), \quad x \in \partial\Omega, t \in [-\tau, 0] \qquad (17)$$

where $p_0^*$ and $s_0^*$ satisfy the correct boundary conditions at $t = 0$ and $s_0(x) = s_-(x, 0)$ on $\partial\Omega$. The definitions of the dimensionless parameters are found in Table 3.

Given the diffusive nature of the system, we anticipate smooth solutions to the system of Equations (11)-(17) that exist globally in time. The proof of existence of such solutions and the properties they enjoy is left as future work. Hereafter, we assume the existence of classical solutions $p, s \in C_x^2 \cap C_t^1(\Omega, [0, \infty))$.

Note that through nondimensionalization and rescaling, we reduce a system with 10 dimensional parameters (Table 1) to a dimensionless system involving 7 dimensionless parameters (Table 3), excluding $m$. It is useful to consider that all dimensional systems resulting in the same set of dimensionless parameters have the same dimensionless solutions.

### 2.1.4 Simplified Geometries

As specific scenarios, we consider $d = 1$ in a planar geometry and $d = 2, 3$ in radially symmetric geometries — effectively describing an idealized “spherical neuron” in different dimensions. See
<table><caption>Table 3: Dimensionless parameters appearing in the model.</caption><thead><tr><th>Parameter</th><th>Definition</th><th>Interpretation</th></tr></thead><tbody><tr><td>γ</td><td>cS<sub>c</sub>/a</td><td>ratio of (the rate P is converted to S when S concentration is S<sub>c</sub>) to (the rate P is cleared)</td></tr><tr><td>σ &ll; 1</td><td>f/a</td><td>ratio of the rate P misfolds to the rate it is cleared</td></tr><tr><td>δ &le; 1</td><td>D<sub>S</sub>/D<sub>P</sub></td><td>ratio of S diffusivity to P diffusivity</td></tr><tr><td>η</td><td>A/S<sub>c</sub>√D<sub>P</sub>a</td><td>ratio of characteristic P-concentration to S<sub>c</sub></td></tr><tr><td>β &le; 1</td><td>b/a</td><td>ratio of rate S is cleared to the rate P is cleared</td></tr><tr><td>τ</td><td>at<sub>d</sub></td><td>delay relative to P-clearance time</td></tr><tr><td>ρ</td><td>R√a/√D<sub>P</sub></td><td>ratio of (the characteristic cell size) to (the characteristic distance P travels before being cleared)</td></tr><tr><td>m &gt; 0</td><td>chosen value</td><td>controls abruptness of switch</td></tr></tbody></table>

with

$$p, s \to 0, \quad |x| \to \infty \qquad (13)$$

$$-\nabla p \cdot \hat{n} = \frac{-1}{1 + (s(\cdot, t - \tau))^m} \quad x \in \partial\Omega, t \ge 0 \qquad (14)$$

$$-\nabla s \cdot \hat{n} = 0, \quad \text{on } \partial\Omega \qquad (15)$$

$$(p(x, 0), s(x, 0)) = (p_{0}^{*}(x), s_{0}^{*}(x)), \quad x \in \bar{\Omega} \qquad (16)$$

$$s(x, t) = s_{-}^{*}(x, t), \quad x \in \partial\Omega, t \in [-\tau, 0] \qquad (17)$$

where $p_0^*$ and $s_0^*$ satisfy the correct boundary conditions at $t = 0$ and $s_0(x) = s_-(x, 0)$ on $\partial\Omega$. The definitions of the dimensionless parameters are found in Table 3.

Given the diffusive nature of the system, we anticipate smooth solutions to the system of Equations (11)-(17) that exist globally in time. The proof of existence of such solutions and the properties they enjoy is left as future work. Hereafter, we assume the existence of classical solutions $p, s \in C_x^2 \cap C_t^1(\Omega, [0, \infty))$.

Note that through nondimensionalization and rescaling, we reduce a system with 10 dimensional parameters (Table 1) to a dimensionless system involving 7 dimensionless parameters (Table 3), excluding $m$. It is useful to consider that all dimensional systems resulting in the same set of dimensionless parameters have the same dimensionless solutions.

### 2.1.4 Simplified Geometries

As specific scenarios, we consider $d = 1$ in a planar geometry and $d = 2, 3$ in radially symmetric geometries — effectively describing an idealized “spherical neuron” in different dimensions. See
Figure 2. We denote $r = |x|$ and reduce $\Omega$ to $B = (\rho, \infty)$. Then

$$ \int_{\Omega} z(x)dx = C_d \int_{\rho}^{\infty} z(r)r^{d-1}dr $$

where

$$ C_d = \begin{cases} 1, & d=1, \\ 2\pi, & d=2, \\ 4\pi, & d=3 \end{cases} $$

is a surface area factor.

In the symmetric geometry, on $B \times (0, \infty)$, we have

$$ p, t = p,rr + \frac{d-1}{r}p, r - \gamma ps - (1 + \sigma)p, \quad (18) $$

$$ s, t = \delta(s,rr + \frac{d-1}{r}s, r) + \eta\gamma ps + \eta\sigma p - \beta s, \quad (19) $$

with

$$ p, s \to 0, \quad r \to \infty, \quad (20) $$

$$ p, r|_{r=\rho} = \frac{-1}{1 + s(\rho, t - \tau)^m} \quad t \ge 0, \quad (21) $$

$$ s, r|_{r=\rho} = 0, \quad (22) $$

$$ (p(r, 0), s(r, 0)) = (p_0(r), s_0(r)), \quad r \in \bar{B}, \quad (23) $$

$$ s(\rho, t) = s_-(t), \quad t \in [-\tau, 0], \quad (24) $$

$$ h_{\text{avg}} = \lim_{\Upsilon \to \infty} \frac{1}{\Upsilon} \int_{0}^{\Upsilon} \left(1 - \frac{1}{1 + s(\rho, t - \tau)^m}\right) dt, \quad (25) $$

$$ \omega_{\text{avg}} = \lim_{\Upsilon \to \infty} \frac{1}{\Upsilon} \times \left| \{t| t \text{ is a strict local maximum for } s(\rho, t - \tau) \text{ on } (0, \Upsilon)\} \right|, \quad (26) $$

where $p_0$ and $s_0$ satisfy the correct boundary conditions at $t = 0$, and where $s_-(0) = s_0(\rho)$. We have added the formulas for the average fraction of time the UPR is activated, Eq. (25), and the average dimensionless frequency of UPR activations, Eq. (26). We shall keep $m = 10$ fixed throughout this work as it allows for a relatively rapid switch term.

We choose initial conditions to model the introduction of a highly localized quantity of S at position $r = r_* > \rho$ so that initially $s \approx 0$ is negligibly small at $r = \rho$. We let

$$ i_p := C_d \int_{\rho}^{\infty} r^{d-1} p_0(r) dr \quad (27) $$

and

$$ i_s := C_d \int_{\rho}^{\infty} r^{d-1} s_0(r) dr \quad (28) $$

represent the initial quantity of the proteins and

$$ -p_{0,r}|_{r=\rho} = n_p \quad (29) $$
Figure 2. We denote $r = |x|$ and reduce $\Omega$ to $B = (\rho, \infty)$. Then

$$\int_{\Omega} z(x)dx = C_d \int_{\rho}^{\infty} z(r)r^{d-1}dr$$

where

$$C_d = \begin{cases} 1, & d=1, \\ 2\pi, & d=2, \\ 4\pi, & d=3 \end{cases}$$

is a surface area factor.

In the symmetric geometry, on $B \times (0, \infty)$, we have

$$p_{,t} = p_{,rr} + \frac{d-1}{r}p_{,r} - \gamma ps - (1 + \sigma)p, \qquad (18)$$

$$s_{,t} = \delta(s_{,rr} + \frac{d-1}{r}s_{,r}) + \eta\gamma ps + \eta\sigma p - \beta s, \qquad (19)$$

with

$$p, s \to 0, \quad r \to \infty, \qquad (20)$$

$$p_{,r}|_{r=\rho} = \frac{-1}{1 + s(\rho, t - \tau)^m} \quad t \ge 0, \qquad (21)$$

$$s_{,r}|_{r=\rho} = 0, \qquad (22)$$

$$(p(r, 0), s(r, 0)) = (p_0(r), s_0(r)), \quad r \in \bar{B}, \qquad (23)$$

$$s(\rho, t) = s_{-}(t), \quad t \in [-\tau, 0], \qquad (24)$$

$$h_{\text{avg}} = \lim_{\Upsilon \to \infty} \frac{1}{\Upsilon} \int_{0}^{\Upsilon} \left(1 - \frac{1}{1 + s(\rho, t - \tau)^m}\right) dt, \qquad (25)$$

$$\omega_{\text{avg}} = \lim_{\Upsilon \to \infty} \frac{1}{\Upsilon} \times |\{t| t \text{ is a strict local maximum for } s(\rho, t - \tau) \text{ on } (0, \Upsilon)\}|, \qquad (26)$$

where $p_0$ and $s_0$ satisfy the correct boundary conditions at $t = 0$, and where $s_-(0) = s_0(\rho)$. We have added the formulas for the average fraction of time the UPR is activated, Eq. (25), and the average dimensionless frequency of UPR activations, Eq. (26). We shall keep $m = 10$ fixed throughout this work as it allows for a relatively rapid switch term.

We choose initial conditions to model the introduction of a highly localized quantity of S at position $r = r_* > \rho$ so that initially $s \approx 0$ is negligibly small at $r = \rho$. We let

$$i_p := C_d \int_{\rho}^{\infty} r^{d-1} p_0(r) dr \qquad (27)$$

and

$$i_s := C_d \int_{\rho}^{\infty} r^{d-1} s_0(r) dr \qquad (28)$$

represent the initial quantity of the proteins and

$$-p_{0,r}|_{r=\rho} = n_p \qquad (29)$$
<figure><img src="image_4.png" /><figcaption>Figure 2: Reductions to 1-dimensional (left), 2-dimensional radially symmetric (middle), and 3-dimensional radially symmetric (right) geometries, respectively. The green dots represent the P proteins, while the red rods represent the S proteins.</figcaption></figure>

specify the normal derivative of $p_0$ at the cell boundary where $s_0 = 0$. We then choose

$$p_0(r) = c_{1d}r^{1-d}e^{-c_{2d}(r-\rho)-(r-\rho)^2}, \qquad (30)$$

$$s_0(r) = \frac{i_s r^{1-d}}{C_d \omega \left(1 + \tanh\left(\frac{r_* - \rho}{\omega}\right)\right)} \mathrm{sech}^2\left(\frac{r - r_*}{\omega}\right), \qquad (31)$$

$$s_-(t) = s_0(\rho), \quad t \in [-\tau, 0], \qquad (32)$$

where $c_{1d}$ and $c_{2d}$ satisfy the system

$$\frac{\sqrt{\pi} c_{1d}}{2} e^{c_{2d}^2/4} \mathrm{erfc}(c_{2d}/2) = i_p/C_d, \\ c_{1d} c_{2d} \rho^{1-d} + (d-1) c_{1d} c_{2d}^d = n_p,$$

and $\omega = 0.2$ is a shape parameter. Under this setup, $C_d \int_{\rho}^{\infty} r^{d-1} s_0(r) dr = i_s$. With $\omega$ small, $s_0$ and its derivatives are very small for $r_* - r = O(1)$ so that $s_0(\rho)$, $s'_0(\rho) \approx 0$ are negligibly small. Then with $s_0(\rho) = 0$, boundary conditions require that $p'_0(\rho) = -1$. The constants $c_{1d}$ and $c_{2d}$ ensure $p'_0(\rho) = -n_p$ and $C_d \int_{\rho}^{\infty} r^{d-1} p_0(r) dr = i_p$. For most problems, $i_p = i_s = n_p = 1$, but this system is more general.

We make one final remark on the total protein quantity. In $d = 1$, we interpret the system as describing protein concentrations in a very narrow, straight region. The cell membrane is one boundary and concentration variations orthogonal to the direction pointing away from the cell are negligible. Thus, the system has a small cross-sectional area and is effectively one-dimensional. Likewise, in $d = 2$, we interpret the system as describing protein concentrations when primarily confined to a plane with limited variation normal to the plane, with concentration varying with distance from a circular cell membrane. The system has a narrow height and is primarily two dimensional with the cell membrane being a circle. In $d = 3$, we interpret the system as concentrations of proteins as their radial distance from the centre of a sphere varies. Thus, for $d = 1$, $i_p$ is expressing a dimensionless total quantity of P per unit area of the region. Likewise with $d = 2$, $i_p$ is the total dimensionless quantity of P per unit height of the region. Then with $d = 3$, $i_p$ is the total dimensionless quantity of P. The same applies to $i_s$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Reductions to 1-dimensional (left), 2-dimensional radially symmetric (middle), and 3-dimensional radially symmetric (right) geometries, respectively. The green dots represent the P proteins, while the red rods represent the S proteins.</figcaption></figure>

specify the normal derivative of $p_0$ at the cell boundary where $s_0 = 0$. We then choose

$$p_0(r) = c_{1d}r^{1-d}e^{-c_{2d}(r-\rho)-(r-\rho)^2}, \qquad (30)$$

$$s_0(r) = \frac{i_s r^{1-d}}{C_d \omega \left(1 + \tanh\left(\frac{r_* - \rho}{\omega}\right)\right)} \mathrm{sech}^2\left(\frac{r - r_*}{\omega}\right), \qquad (31)$$

$$s_-(t) = s_0(\rho), \quad t \in [-\tau, 0], \qquad (32)$$

where $c_{1d}$ and $c_{2d}$ satisfy the system

$$\frac{\sqrt{\pi} c_{1d}}{2} e^{c_{2d}^2/4} \mathrm{erfc}(c_{2d}/2) = i_p/C_d, \\ c_{1d} c_{2d} \rho^{1-d} + (d-1) c_{1d} c_{2d}^d = n_p,$$

and $\omega = 0.2$ is a shape parameter. Under this setup, $C_d \int_{\rho}^{\infty} r^{d-1} s_0(r) dr = i_s$. With $\omega$ small, $s_0$ and its derivatives are very small for $r_* - r = O(1)$ so that $s_0(\rho)$, $s'_0(\rho) \approx 0$ are negligibly small. Then with $s_0(\rho) = 0$, boundary conditions require that $p'_0(\rho) = -1$. The constants $c_{1d}$ and $c_{2d}$ ensure $p'_0(\rho) = -n_p$ and $C_d \int_{\rho}^{\infty} r^{d-1} p_0(r) dr = i_p$. For most problems, $i_p = i_s = n_p = 1$, but this system is more general.

We make one final remark on the total protein quantity. In $d = 1$, we interpret the system as describing protein concentrations in a very narrow, straight region. The cell membrane is one boundary and concentration variations orthogonal to the direction pointing away from the cell are negligible. Thus, the system has a small cross-sectional area and is effectively one-dimensional. Likewise, in $d = 2$, we interpret the system as describing protein concentrations when primarily confined to a plane with limited variation normal to the plane, with concentration varying with distance from a circular cell membrane. The system has a narrow height and is primarily two dimensional with the cell membrane being a circle. In $d = 3$, we interpret the system as concentrations of proteins as their radial distance from the centre of a sphere varies. Thus, for $d = 1$, $i_p$ is expressing a dimensionless total quantity of P per unit area of the region. Likewise with $d = 2$, $i_p$ is the total dimensionless quantity of P per unit height of the region. Then with $d = 3$, $i_p$ is the total dimensionless quantity of P. The same applies to $i_s$.
## 2.2 Numerical Method

Owing to the challenging nature of the equations, we analyze our model numerically. Full details can be found in the Appendix A. In brief: we implement a finite difference scheme that is first order in time, second order in space, with semi-implicit timestepping.

# 3 Math and Model Investigations

This section is dedicated to studying the qualitative behavior of the model itself over a range of parameters. Numerical specifications are given in Appendix B.1.

## 3.1 Model Simulations

To study how the system responds to different parameters, we begin by considering variations in $\beta$ and $\delta$ at different fixed values of $\tau$. In making these explorations, we focus primarily on phase portraits for $d = 1$ dimensions depicting the protein concentrations at the cell membrane over time (since this is the driving mechanism of the UPR). This is found in Figures 3-7. We also consider the spatiotemporal spreading of the proteins with further experiments in Figures 8 and 9. We wish to caution the reader that the ranges of values along the vertical differ between trials.

From the phase portraits, Figures 3-7, we observe that oscillations tend to decrease in intensity and are dampened out as $\beta$ increases, as $\delta$ increases, and as $\tau$ decreases. The UPR intensity $h_{\text{avg}}$ tends to decrease as $\beta$ and $\delta$ increase, and does not vary significantly as $\tau$ varies in the experiments. The UPR frequency $\omega_{\text{avg}}$, when present, is quite consistent in its values across the phase portraits for fixed $\tau$-values, and tends to decrease as $\tau$ increases.

From the spatiotemporal plots for parameter variations, Figures 8-9, we observe that the oscillations become less pronounced as the dimension $d$ increases. We again observe how decreasing $\delta$ and $\beta$ can make the system oscillatory. Increasing $\gamma$ makes these oscillations more pronounced. The effect of $\sigma$ on the oscillations is quite small. Increasing $\eta$ makes the oscillations more pronounced in the $d=2$ case and drastically increases the $s$-values in $d=3$. Finally, when $\rho$ is increased, the $d=2$ system becomes more oscillatory and the $d=3$ system begins to oscillate (the only set of oscillations observed in the spatiotemporal plots). The fact that $\rho$ has no effect in $d=1$ is unsurprising since there is no geometric $\frac{d-1}{r}$ term appearing in the Laplacian. As $\rho$ increases, the $d=2$ and $d=3$ systems are more planar on the cell membrane, so it makes sense that with $d=1$ being highly oscillatory in general, increasing $\rho$ should make $d=2$ and $d=3$ more oscillatory.

In Figure 10, we illustrate how the oscillations and UPR intensity are weakened as the dimension of the system $d$ increases.

## 3.2 Comparison with Previous Work

Having explored many of the parameters in our model, we briefly turn to the question of how this new model compares with the DDE model of Adimy [22] which inspired this work. The DDE model also studied a heterodimer system with a production rate for healthy monomers with delay, clearance rates for both the healthy and toxic proteins, and a recruitment rate to convert healthy proteins to toxic proteins. The authors found that for low levels of delay, the system could be stable (without oscillations) and that a bifurcation could take place by increasing the delay to yield an oscillatory solution. With the DDEs, a locally and globally asymptotically stable disease-free steady
## 2.2 Numerical Method

Owing to the challenging nature of the equations, we analyze our model numerically. Full details can be found in the Appendix A. In brief: we implement a finite difference scheme that is first order in time, second order in space, with semi-implicit timestepping.

# 3 Math and Model Investigations

This section is dedicated to studying the qualitative behavior of the model itself over a range of parameters. Numerical specifications are given in Appendix B.1.

## 3.1 Model Simulations

To study how the system responds to different parameters, we begin by considering variations in $\beta$ and $\delta$ at different fixed values of $\tau$. In making these explorations, we focus primarily on phase portraits for $d = 1$ dimensions depicting the protein concentrations at the cell membrane over time (since this is the driving mechanism of the UPR). This is found in Figures 3-7. We also consider the spatiotemporal spreading of the proteins with further experiments in Figures 8 and 9. We wish to caution the reader that the ranges of values along the vertical differ between trials.

From the phase portraits, Figures 3-7, we observe that oscillations tend to decrease in intensity and are dampened out as $\beta$ increases, as $\delta$ increases, and as $\tau$ decreases. The UPR intensity $h_{\text{avg}}$ tends to decrease as $\beta$ and $\delta$ increase, and does not vary significantly as $\tau$ varies in the experiments. The UPR frequency $\omega_{\text{avg}}$, when present, is quite consistent in its values across the phase portraits for fixed $\tau$-values, and tends to decrease as $\tau$ increases.

From the spatiotemporal plots for parameter variations, Figures 8-9, we observe that the oscillations become less pronounced as the dimension $d$ increases. We again observe how decreasing $\delta$ and $\beta$ can make the system oscillatory. Increasing $\gamma$ makes these oscillations more pronounced. The effect of $\sigma$ on the oscillations is quite small. Increasing $\eta$ makes the oscillations more pronounced in the $d=2$ case and drastically increases the $s$-values in $d=3$. Finally, when $\rho$ is increased, the $d=2$ system becomes more oscillatory and the $d=3$ system begins to oscillate (the only set of oscillations observed in the spatiotemporal plots). The fact that $\rho$ has no effect in $d=1$ is unsurprising since there is no geometric $\frac{d-1}{r}$ term appearing in the Laplacian. As $\rho$ increases, the $d=2$ and $d=3$ systems are more planar on the cell membrane, so it makes sense that with $d=1$ being highly oscillatory in general, increasing $\rho$ should make $d=2$ and $d=3$ more oscillatory.

In Figure 10, we illustrate how the oscillations and UPR intensity are weakened as the dimension of the system $d$ increases.

## 3.2 Comparison with Previous Work

Having explored many of the parameters in our model, we briefly turn to the question of how this new model compares with the DDE model of Adimy [22] which inspired this work. The DDE model also studied a heterodimer system with a production rate for healthy monomers with delay, clearance rates for both the healthy and toxic proteins, and a recruitment rate to convert healthy proteins to toxic proteins. The authors found that for low levels of delay, the system could be stable (without oscillations) and that a bifurcation could take place by increasing the delay to yield an oscillatory solution. With the DDEs, a locally and globally asymptotically stable disease-free steady
<figure><img src="image_5.png" /><figcaption>Figure 3: Phase diagram at $\tau = 0$ with $d = 1$ along membrane $r = \rho$. The horizontal axis is $t$, the vertical axis is concentration with $p$ solid blue and $s$ dashed red. Unspecified parameter values are $\gamma = 1$, $\sigma = 0.02$, $\eta = 1$, and $\rho = 0.25$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Phase diagram at $\tau = 0$ with $d = 1$ along membrane $r = \rho$. The horizontal axis is $t$, the vertical axis is concentration with $p$ solid blue and $s$ dashed red. Unspecified parameter values are $\gamma = 1$, $\sigma = 0.02$, $\eta = 1$, and $\rho = 0.25$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Phase diagram at $\tau = 1$ with $d = 1$ along membrane $r = \rho$. The horizontal axis is $t$, the vertical axis is concentration with $p$ solid blue and $s$ dashed red. Unspecified parameter values are $\gamma = 1$, $\sigma = 0.02$, $\eta = 1$, and $\rho = 0.25$.</figcaption></figure>
<table><caption>Figure 4: Phase diagram at $τ = 1$ with $d = 1$ along membrane $r = ρ$. The horizontal axis is $t$, the vertical axis is concentration with $p$ solid blue and $s$ dashed red. Unspecified parameter values are $γ = 1$, $σ = 0.02$, $η = 1$, and $ρ = 0.25$.</caption><thead><tr><th>δ</th><th>0.05</th><th>0.1</th><th>0.15</th><th>0.2</th><th>0.25</th><th>0.3</th></tr></thead><tbody><tr><th>0.05</th><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>0.1</th><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>0.15</th><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>0.2</th><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>0.25</th><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>0.3</th><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>Figure 5: Phase diagram at $\tau = 2$ with $d = 1$ along membrane $r = \rho$. The horizontal axis is $t$, the vertical axis is concentration with $p$ solid blue and $s$ dashed red. Unspecified parameter values are $\gamma = 1$, $\sigma = 0.02$, $\eta = 1$, and $\rho = 0.25$.</caption><thead><tr><th>δ</th><th>0.05</th><th>0.1</th><th>0.15</th><th>0.2</th><th>0.25</th><th>0.3</th></tr></thead><tbody><tr><th>0.05</th><td>2.2 β = 0.758 ω<sub>avg</sub> = 0.123</td><td>2.1 β = 0.721 ω<sub>avg</sub> = 0.125</td><td>1.8 β = 0.689 ω<sub>avg</sub> = 0.127</td><td>1.7 β = 0.666 ω<sub>avg</sub> = 0.128</td><td>1.6 β = 0.644 ω<sub>avg</sub> = 0.128</td><td>1.6 β = 0.623 ω<sub>avg</sub> = 0.129</td></tr><tr><th>0.1</th><td>2 β = 0.643 ω<sub>avg</sub> = 0.12</td><td>1.7 β = 0.58 ω<sub>avg</sub> = 0.124</td><td>1.6 β = 0.549 ω<sub>avg</sub> = 0.126</td><td>1.6 β = 0.522 ω<sub>avg</sub> = 0.126</td><td>1.6 β = 0.495 ω<sub>avg</sub> = 0.129</td><td>1.6 β = 0.468 ω<sub>avg</sub> = 0.129</td></tr><tr><th>0.15</th><td>1.8 β = 0.54 ω<sub>avg</sub> = 0.122</td><td>1.6 β = 0.476 ω<sub>avg</sub> = 0.126</td><td>1.6 β = 0.436 ω<sub>avg</sub> = 0.127</td><td>1.6 β = 0.4 ω<sub>avg</sub> = 0.127</td><td>1.6 β = 0.365 ω<sub>avg</sub> = 0.128</td><td>1.6 β = 0.332 ω<sub>avg</sub> = 0.128</td></tr><tr><th>0.2</th><td>1.6 β = 0.452 ω<sub>avg</sub> = 0.123</td><td>1.6 β = 0.38 ω<sub>avg</sub> = 0.124</td><td>1.6 β = 0.327 ω<sub>avg</sub> = 0.125</td><td>1.6 β = 0.282 ω<sub>avg</sub> = 0.126</td><td>1.6 β = 0.241 ω<sub>avg</sub> = 0.127</td><td>1.6 β = 0.205 ω<sub>avg</sub> = 0.126</td></tr><tr><th>0.25</th><td>1.6 β = 0.382 ω<sub>avg</sub> = 0.122</td><td>1.6 β = 0.298 ω<sub>avg</sub> = 0.122</td><td>1.6 β = 0.229 ω<sub>avg</sub> = 0.124</td><td>1.6 β = 0.177 ω<sub>avg</sub> = 0.124</td><td>1.6 β = 0.136 ω<sub>avg</sub> = 0.119</td><td>1.6 β = 0.102 ω<sub>avg</sub> = 0.111</td></tr><tr><th>0.3</th><td>1.6 β = 0.306 ω<sub>avg</sub> = 0.119</td><td>1.6 β = 0.207 ω<sub>avg</sub> = 0.123</td><td>1.6 β = 0.136 ω<sub>avg</sub> = 0.121</td><td>1.6 β = 0.089 ω<sub>avg</sub> = 0.11</td><td>1.6 β = 0.053 ω<sub>avg</sub> = 0</td><td>1.6 β = 0.028 ω<sub>avg</sub> = 0</td></tr></tbody></table>
<table><caption>Figure 5: Phase diagram at $\tau = 2$ with $d = 1$ along membrane $r = \rho$. The horizontal axis is $t$, the vertical axis is concentration with $p$ solid blue and $s$ dashed red. Unspecified parameter values are $\gamma = 1$, $\sigma = 0.02$, $\eta = 1$, and $\rho = 0.25$.</caption><thead><tr><th>δ</th><th>0.05</th><th>0.1</th><th>0.15</th><th>0.2</th><th>0.25</th><th>0.3</th></tr></thead><tbody><tr><th>0.05</th><td>2.2 β = 0.758 ω<sub>avg</sub> = 0.123</td><td>2.2 β = 0.721 ω<sub>avg</sub> = 0.125</td><td>1.8 β = 0.689 ω<sub>avg</sub> = 0.127</td><td>1.7 β = 0.666 ω<sub>avg</sub> = 0.128</td><td>1.6 β = 0.644 ω<sub>avg</sub> = 0.128</td><td>1.6 β = 0.623 ω<sub>avg</sub> = 0.129</td></tr><tr><th>0.1</th><td>2 β = 0.643 ω<sub>avg</sub> = 0.12</td><td>1.7 β = 0.58 ω<sub>avg</sub> = 0.124</td><td>1.6 β = 0.549 ω<sub>avg</sub> = 0.126</td><td>1.6 β = 0.522 ω<sub>avg</sub> = 0.126</td><td>1.6 β = 0.495 ω<sub>avg</sub> = 0.129</td><td>1.6 β = 0.468 ω<sub>avg</sub> = 0.129</td></tr><tr><th>0.15</th><td>1.8 β = 0.54 ω<sub>avg</sub> = 0.122</td><td>1.6 β = 0.476 ω<sub>avg</sub> = 0.126</td><td>1.6 β = 0.436 ω<sub>avg</sub> = 0.127</td><td>1.6 β = 0.4 ω<sub>avg</sub> = 0.127</td><td>1.6 β = 0.365 ω<sub>avg</sub> = 0.128</td><td>1.6 β = 0.332 ω<sub>avg</sub> = 0.128</td></tr><tr><th>0.2</th><td>1.6 β = 0.452 ω<sub>avg</sub> = 0.123</td><td>1.6 β = 0.38 ω<sub>avg</sub> = 0.124</td><td>1.6 β = 0.327 ω<sub>avg</sub> = 0.125</td><td>1.6 β = 0.282 ω<sub>avg</sub> = 0.126</td><td>1.6 β = 0.241 ω<sub>avg</sub> = 0.127</td><td>1.6 β = 0.205 ω<sub>avg</sub> = 0.126</td></tr><tr><th>0.25</th><td>1.6 β = 0.382 ω<sub>avg</sub> = 0.122</td><td>1.6 β = 0.298 ω<sub>avg</sub> = 0.122</td><td>1.6 β = 0.229 ω<sub>avg</sub> = 0.124</td><td>1.6 β = 0.177 ω<sub>avg</sub> = 0.124</td><td>1.6 β = 0.136 ω<sub>avg</sub> = 0.119</td><td>1.6 β = 0.102 ω<sub>avg</sub> = 0.113</td></tr><tr><th>0.3</th><td>1.6 β = 0.306 ω<sub>avg</sub> = 0.119</td><td>1.6 β = 0.207 ω<sub>avg</sub> = 0.123</td><td>1.6 β = 0.136 ω<sub>avg</sub> = 0.121</td><td>1.6 β = 0.089 ω<sub>avg</sub> = 0.11</td><td>1.6 β = 0.053 ω<sub>avg</sub> = 0</td><td>1.6 β = 0.028 ω<sub>avg</sub> = 0</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 6: Phase diagram at $\tau = 3$ with $d = 1$ along membrane $r = \rho$. The horizontal axis is $t$, the vertical axis is concentration with $p$ solid blue and $s$ dashed red. Unspecified parameter values are $\gamma = 1$, $\sigma = 0.02$, $\eta = 1$, and $\rho = 0.25$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 6: Phase diagram at $\tau = 3$ with $d = 1$ along membrane $r = \rho$. The horizontal axis is $t$, the vertical axis is concentration with $p$ solid blue and $s$ dashed red. Unspecified parameter values are $\gamma = 1$, $\sigma = 0.02$, $\eta = 1$, and $\rho = 0.25$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Phase diagram at $\tau = 4$ with $d = 1$ along membrane $r = \rho$. The horizontal axis is $t$, the vertical axis is concentration with $p$ solid blue and $s$ dashed red. Unspecified parameter values are $\gamma = 1$, $\sigma = 0.02$, $\eta = 1$, and $\rho = 0.25$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 7: Phase diagram at $\tau = 4$ with $d = 1$ along membrane $r = \rho$. The horizontal axis is $t$, the vertical axis is concentration with $p$ solid blue and $s$ dashed red. Unspecified parameter values are $\gamma = 1$, $\sigma = 0.02$, $\eta = 1$, and $\rho = 0.25$.</figcaption></figure>
<table><thead><tr><th>p Behavior</th><th>Dimension 1</th><th>Dimension 2</th><th>Dimension 3</th></tr></thead><tbody><tr><td>γ = 1<br>σ = 0<br>δ = β = 0.3<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 1<br>σ = 0<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 1<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 1<br>τ = 2.5<br>ρ = 0.5</td><td></td><td></td><td></td></tr></tbody></table>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

<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" /><figcaption>Figure 8: Behavior of *p* with varying parameters.</figcaption></figure>
<table><thead><tr><th>p Behavior</th><th>Dimension 1</th><th>Dimension 2</th><th>Dimension 3</th></tr></thead><tbody><tr><td>γ = 1<br>σ = 0<br>δ = β = 0.3<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 1<br>σ = 0<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 1<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 1<br>τ = 2.5<br>ρ = 0.5</td><td></td><td></td><td></td></tr></tbody></table>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

<figure><img src="image_24.png" /><figcaption>Figure 8: Behavior of *p* with varying parameters.</figcaption></figure>
<table><thead><tr><th>s Behavior</th><th>Dimension 1</th><th>Dimension 2</th><th>Dimension 3</th></tr></thead><tbody><tr><td>γ = 1<br>σ = 0<br>δ = β = 0.3<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 1<br>σ = 0<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 1<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 1<br>τ = 2.5<br>ρ = 0.5</td><td></td><td></td><td></td></tr></tbody></table>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

<figure><img src="image_44.png" /><figcaption>Figure 9: Behavior of s with varying parameters.</figcaption></figure>
<table><thead><tr><th>s Behavior</th><th>Dimension 1</th><th>Dimension 2</th><th>Dimension 3</th></tr></thead><tbody><tr><td>γ = 1<br>σ = 0<br>δ = β = 0.3<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 1<br>σ = 0<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 0.5<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 1<br>τ = 2.5<br>ρ = 0.25</td><td></td><td></td><td></td></tr><tr><td>γ = 5<br>σ = 0.05<br>δ = β = 0.05<br>η = 1<br>τ = 2.5<br>ρ = 0.5</td><td></td><td></td><td></td></tr></tbody></table>

<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" /></figure>

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

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

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

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

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

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

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

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

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

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

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

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

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

<figure><img src="image_42.png" /><figcaption>Figure 9: Behavior of s with varying parameters.</figcaption></figure>
<figure><img src="image_45.png" /><figcaption>Figure 10: Temporal variations in concentration at $r = \rho$ in $d = 1$ (left), $d = 2$ (middle), and $d = 3$ (right) with $p$ solid blue and $s$ dashed red. Parameters are $\gamma = 6$, $\sigma = 0$, $\delta = \beta = 0.075$, $\eta = 0.9$, $\tau = 1.25$, and $\rho = 0.5$.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>Figure 11: Temporal variations in concentration at $r = \rho$ with $d = 1$ with $p$ solid blue and $s$ dashed red. From left to right: disease-free ($\gamma = 0.2, \sigma = 0, \delta = 0.2, \eta = 0.5, \beta = 1, \tau = 0.2, \text{and } \rho = 0.25$); endemic, no oscillations ($\gamma = 3, \sigma = 0, \delta = 0.2, \eta = 0.5, \beta = 0.2, \tau = 0.2, \text{and } \rho = 0.25$); endemic, oscillations ($\gamma = 3, \sigma = 0, \delta = 0.2, \eta = 0.5, \beta = 0.2, \tau = 2, \text{and } \rho = 0.25$); and endemic, no oscillations ($\gamma = 3, \sigma = 0, \delta = 1, \eta = 0.5, \beta = 0.2, \tau = 2, \text{and } \rho = 0.25$).</figcaption></figure>

state was found when the reproductive ratio $R_0 < 1$, where $R_0$ gives the ratio of the product of the peak production rate and recruitment rate to the product of the clearance rates. When $R_0 > 1$, the disease takes over and oscillatory solutions may occur.

Here, we illustrate how our system observes similar behavior. In Figure 11, we begin with parameters for a disease-free state with small $\gamma$ and large $\beta$, bring about a disease state with increasing $\gamma$, induce oscillations with an increase in $\tau$, and then eliminate the oscillations with an increase in $\delta$. We note that while the behavior is similar to the DDE model, factors like $\delta$ stemming from diffusion also play a role. And as noted in Figures 8-10, geometric factors such as the cell size $\rho$ and the dimension $d$ of the spatial domain are relevant.

## 3.3 Mathematical Perspective

Given the nonlinear structure in the model, especially within the delayed flux boundary condition for $P$, there are open questions as to the well-posedness of the PDE system and to the solution properties. For example, a function in $g \in C_x^2 C_t^1(\Omega; (0, \infty))$ could exhibit unbounded behavior for $x$ near $\partial\Omega$ as $t \to \infty$ while retaining a bounded normal derivative. We expect solutions with intuitive, physical bounds to exist, but this is left as a future work in mathematics. Additionally, while delay differential equations are well-studied, having a delay in the boundary condition itself is somewhat unique.
<figure><img src="image_43.png" /><figcaption>Figure 10: Temporal variations in concentration at $r = \rho$ in $d = 1$ (left), $d = 2$ (middle), and $d = 3$ (right) with $p$ solid blue and $s$ dashed red. Parameters are $\gamma = 6$, $\sigma = 0$, $\delta = \beta = 0.075$, $\eta = 0.9$, $\tau = 1.25$, and $\rho = 0.5$.</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Figure 11: Temporal variations in concentration at $r = \rho$ with $d = 1$ with $p$ solid blue and $s$ dashed red. From left to right: disease-free ($\gamma = 0.2, \sigma = 0, \delta = 0.2, \eta = 0.5, \beta = 1, \tau = 0.2, \text{and } \rho = 0.25$); endemic, no oscillations ($\gamma = 3, \sigma = 0, \delta = 0.2, \eta = 0.5, \beta = 0.2, \tau = 0.2, \text{and } \rho = 0.25$); endemic, oscillations ($\gamma = 3, \sigma = 0, \delta = 0.2, \eta = 0.5, \beta = 0.2, \tau = 2, \text{and } \rho = 0.25$); and endemic, no oscillations ($\gamma = 3, \sigma = 0, \delta = 1, \eta = 0.5, \beta = 0.2, \tau = 2, \text{and } \rho = 0.25$).</figcaption></figure>

state was found when the reproductive ratio $R_0 < 1$, where $R_0$ gives the ratio of the product of the peak production rate and recruitment rate to the product of the clearance rates. When $R_0 > 1$, the disease takes over and oscillatory solutions may occur.

Here, we illustrate how our system observes similar behavior. In Figure 11, we begin with parameters for a disease-free state with small $\gamma$ and large $\beta$, bring about a disease state with increasing $\gamma$, induce oscillations with an increase in $\tau$, and then eliminate the oscillations with an increase in $\delta$. We note that while the behavior is similar to the DDE model, factors like $\delta$ stemming from diffusion also play a role. And as noted in Figures 8-10, geometric factors such as the cell size $\rho$ and the dimension $d$ of the spatial domain are relevant.

## 3.3 Mathematical Perspective

Given the nonlinear structure in the model, especially within the delayed flux boundary condition for $P$, there are open questions as to the well-posedness of the PDE system and to the solution properties. For example, a function in $g \in C_x^2 C_t^1(\Omega; (0, \infty))$ could exhibit unbounded behavior for $x$ near $\partial\Omega$ as $t \to \infty$ while retaining a bounded normal derivative. We expect solutions with intuitive, physical bounds to exist, but this is left as a future work in mathematics. Additionally, while delay differential equations are well-studied, having a delay in the boundary condition itself is somewhat unique.
We additionally observe several “bifurcations” that emerge in our system, such as the turning off/on of oscillations as model parameters vary or going from a disease-free to diseased state. We have 7 dimensionless parameters (excluding *m*) and it would be of great interest to understand their individual effects from a rigorous, theoretical perspective. We also note that, due to the limitations of a numerical study, some of the observed oscillations that appear to be damped may actually be sustained at a small amplitude in the limit, and some of the systems that appear to have no oscillations may actually oscillate. This makes further mathematical study important.

# 4 Biological Implications

We now focus our attention on the biological and medical insights provided by the model.

## 4.1 Parameter Estimation and Setup

To our knowledge, extant literature does not provide data to accurately determine the many parameters used in our model. It would be of interest to conduct experiments that directly measure the chemical reaction rates, diffusivities, etc., in vivo. Nevertheless, we venture to obtain very loose estimates for the parameters based on the literature to then model, as best possible, a representative biological system. The data we use come from a variety of animal models, tissue types, and mathematical models with different objectives than our own. Our parameters are listed in Table 4 in cgs units and justifications follow in the paragraphs below. The resulting dimensionless parameters are then given in Table 5. We only utilize one significant figure to reflect the parameter uncertainty.

**Parameters $D_P$ and $D_S$:** The cellular prion protein $PrP^C$ has a molecular mass of 27–30 kDa [33]. We shall use 28.5 kDa. Interstitial fluid has been reported to have a range of viscosities of 0.7–3.5 cP [34] and we shall use the mean of the extremes for a value 2.1 cP. A useful engineering approximation is that the diffusivity of a protein is given by

$$D \approx \frac{\zeta T}{vm^{1/3}} \qquad (33)$$

where $\zeta = 8.34 \times 10^{-8} \text{ cm}^2 \cdot \text{cP} \cdot \text{Da}^{1/3} \text{ s} \cdot (\circ\text{K})$, $T$ is the absolute temperature, $v$ is the solvent viscosity, and $m$ is the molar mass of the solute [35]. However, due to tortuosity, this is scaled down further by $\iota^2$ where $\iota \approx 1.6$ is the tortuosity of the brain [36]. At $T = 310^\circ\text{K}$, this yields $D_P \approx 2 \times 10^{-7} \text{ cm}^2/\text{s}$.

In a nucleated polymerization model based on mouse brain data, the mean number of monomers in a $PrP^Sc$ protein is estimated at 100–10000 [37]. Using the geometric mean, this suggests our S-proteins should have a mass that is 1000 times larger than our P-proteins. Based on diffusivity scaling with the inverse cube root of molecular mass in Equation (33), we find $D_S \approx 2 \times 10^{-8} \text{ cm}^2/\text{s}$.

**Parameters $a$ and $b$:** Based on mouse and hamster brain data, a nucleated polymerization model estimated the clearance rate of $PrP^C$ at approximately 3–5 day$^{-1}$ [37]. We use a value of 4 day$^{-1}$ so $a \approx 5 \times 10^{-5} \text{ s}^{-1}$. The rate of clearance of toxic Scrapie proteins was also estimated to be 0.03–0.2 day$^{-1}$ and we use the geometric mean of the extremes yielding $b \approx 0.0775 \text{ day}^{-1} \approx 9 \times 10^{-7} \text{ s}^{-1}$.
We additionally observe several “bifurcations” that emerge in our system, such as the turning off/on of oscillations as model parameters vary or going from a disease-free to diseased state. We have 7 dimensionless parameters (excluding *m*) and it would be of great interest to understand their individual effects from a rigorous, theoretical perspective. We also note that, due to the limitations of a numerical study, some of the observed oscillations that appear to be damped may actually be sustained at a small amplitude in the limit, and some of the systems that appear to have no oscillations may actually oscillate. This makes further mathematical study important.

# 4 Biological Implications

We now focus our attention on the biological and medical insights provided by the model.

## 4.1 Parameter Estimation and Setup

To our knowledge, extant literature does not provide data to accurately determine the many parameters used in our model. It would be of interest to conduct experiments that directly measure the chemical reaction rates, diffusivities, etc., in vivo. Nevertheless, we venture to obtain very loose estimates for the parameters based on the literature to then model, as best possible, a representative biological system. The data we use come from a variety of animal models, tissue types, and mathematical models with different objectives than our own. Our parameters are listed in Table 4 in cgs units and justifications follow in the paragraphs below. The resulting dimensionless parameters are then given in Table 5. We only utilize one significant figure to reflect the parameter uncertainty.

**Parameters $D_P$ and $D_S$:** The cellular prion protein $PrP^C$ has a molecular mass of 27–30 kDa [33]. We shall use 28.5 kDa. Interstitial fluid has been reported to have a range of viscosities of 0.7–3.5 cP [34] and we shall use the mean of the extremes for a value 2.1 cP. A useful engineering approximation is that the diffusivity of a protein is given by

$$D \approx \frac{\zeta T}{vm^{1/3}} \qquad (33)$$

where $\zeta = 8.34 \times 10^{-8} \text{ cm}^2 \cdot \text{cP} \cdot \text{Da}^{1/3} \text{ s} \cdot (\circ \text{K})$, $T$ is the absolute temperature, $v$ is the solvent viscosity, and $m$ is the molar mass of the solute [35]. However, due to tortuosity, this is scaled down further by $\iota^2$ where $\iota \approx 1.6$ is the tortuosity of the brain [36]. At $T = 310^\circ \text{K}$, this yields $D_P \approx 2 \times 10^{-7} \text{ cm}^2/\text{s}$.

In a nucleated polymerization model based on mouse brain data, the mean number of monomers in a $PrP^Sc$ protein is estimated at 100–10000 [37]. Using the geometric mean, this suggests our S-proteins should have a mass that is 1000 times larger than our P-proteins. Based on diffusivity scaling with the inverse cube root of molecular mass in Equation (33), we find $D_S \approx 2 \times 10^{-8} \text{ cm}^2/\text{s}$.

**Parameters $a$ and $b$:** Based on mouse and hamster brain data, a nucleated polymerization model estimated the clearance rate of $PrP^C$ at approximately 3–5 day$^{-1}$ [37]. We use a value of 4 day$^{-1}$ so $a \approx 5 \times 10^{-5} \text{ s}^{-1}$. The rate of clearance of toxic Scrapie proteins was also estimated to be 0.03–0.2 day$^{-1}$ and we use the geometric mean of the extremes yielding $b \approx 0.0775 \text{ day}^{-1} \approx 9 \times 10^{-7} \text{ s}^{-1}$.
**Parameter S<sub>c</sub>:** For hamsters at the terminal stage of Scrapie, the concentration of Scrapie proteins was measured to be about 100 µg per gram of brain tissue [38]. The density of brain tissue is approximately 1 g/cm<sup>3</sup> [39]. Taking the terminal concentration as the critical concentration S<sub>c</sub>, we have S<sub>c</sub> = 100 µg/cm<sup>3</sup> = 1 × 10<sup>-4</sup> g/cm<sup>3</sup>.

**Parameter f:** Due to the rarity of spontaneous prion diseases, we take f = 0 s<sup>-1</sup>.

**Parameter c:** Based on hamster data, a heterodimer model estimates the combination rate at 0.15 $\frac{g_{\text{brain}}}{\mu\text{g}\cdot\text{day}}$ [40]. Again with a brain density of 1 g/cm<sup>3</sup> [39], we find c ≈ 0.15 $\frac{\text{cm}^3}{\mu\text{g}\cdot\text{day}}$ ≈ 2$\frac{\text{cm}^3}{\text{g}\cdot\text{s}}$.

**Parameter t<sub>d</sub>:** A study that infected human fibrosarcoma cells with Dengue Virus found an increase in phosphorylated eIF2α — indicating the activation of Integrated Stress Response arm of the UPR — at approximately 6 hours post infection [41]. At 24 – 36 hours in that same study, the inositol-requiring-protein-1 UPR pathway was activated. This suggests a timescale for the UPR may be in the loose range of 6 – 36 hours, and we use the mean of the extremes, t<sub>d</sub> ≈ 21 hours or t<sub>d</sub> ≈ 80000 s.

**Parameter A:** The most challenging parameter to estimate is the maximum flux of P-proteins. We first note that a heterodimer model on hamster data estimates the rate of monomer production to be $\lambda = 4 \frac{\mu g}{(g_{\text{brain}})\cdot\text{day}}$ [40], which with brain density of 1 g/cm<sup>3</sup> [39], we have $\lambda \approx 5 \times 10^{-11} \frac{g}{\text{cm}^3\text{s}}$.

There are approximately 2 × 10<sup>9</sup> neurons per cm<sup>3</sup> [42] across many species. Neuronal perikaryal diameters range from about 6 to 80 µm and the perikarya accounts for approximately 10% of the neuron surface area [43]. By these estimate and with the geometric mean of the diameter lower and upper bounds taken as a length scale, the surface area of a neuron is ~ 5000 µm<sup>2</sup> or 5 × 10<sup>-5</sup> cm<sup>2</sup>.

Then, as a ballpark figure for the flux of P-proteins from a single neuron, we obtain $A^* \approx \frac{5 \times 10^{-11}}{2 \times 10^9 \times 5 \times 10^{-5}} = 5 \times 10^{-16} \frac{g}{\text{cm}^2 \cdot \text{s}}$. To this figure, we make a scaling correction. As our intent with the model is to understand the dynamics experienced by a single typical neuron, we need to take into account that no neuron functions in isolation and the brain is full of neurons. Neurons in a nearby vicinity will also experience similar stresses due to the presence of S-proteins. We thereby choose to scale $A^*$ by a ballpark estimate of the number of neurons that could influence any particular neuron through their release of P-proteins.

On average, a P-protein survives a time of $1/a$, resulting in a characteristic diffusive lengthscale of $\sqrt{D_p}/a \approx 0.06$ cm. Within a volume of 0.063 cm<sup>3</sup>, there are approximately 4 × 10<sup>5</sup> neurons, and we scale $A^*$ up by this to give a nominal flux of $A \approx 2 \times 10^{-10} \frac{g}{\text{cm}^2 \cdot \text{s}}$.

**Parameter R:** cells of the globus pallidus in humans can have dendritic extends of ≈ 1000 µm and soma radii of ≈ 16.5 µm [44]. To come up with a representative length scale for an idealized “spherical neuron”, we take a geometric mean of these length scales giving R ~ 128 µm ≈ 0.01 cm.

**Setup for Systems** We focus on the behavior of our model with the estimated biological parameters. Numerical specifications for our biologically-driven study are given in Appendix B.2. The parameters of Table 4 define characteristic scales of $\bar{t} = 20000$ s, $\bar{x} = 0.06$ cm, $\bar{P} = 6 \times 10^{-5}$ g/cm<sup>3</sup>, and $\bar{S} = 0.0001$ g/cm<sup>3</sup>. Throughout all simulations, we choose initial conditions so that the
**Parameter S<sub>c</sub>:** For hamsters at the terminal stage of Scrapie, the concentration of Scrapie proteins was measured to be about 100 µg per gram of brain tissue [38]. The density of brain tissue is approximately 1 g/cm<sup>3</sup> [39]. Taking the terminal concentration as the critical concentration S<sub>c</sub>, we have S<sub>c</sub> = 100 µg/cm<sup>3</sup> = 1 × 10<sup>-4</sup> g/cm<sup>3</sup>.

**Parameter f:** Due to the rarity of spontaneous prion diseases, we take $f = 0 \text{ s}^{-1}$.

**Parameter c:** Based on hamster data, a heterodimer model estimates the combination rate at 0.15 $\frac{\text{g brain}}{\mu\text{g}\cdot\text{day}}$ [40]. Again with a brain density of 1 g/cm<sup>3</sup> [39], we find $c \approx 0.15 \frac{\text{cm}^3}{\mu\text{g}\cdot\text{day}} \approx 2 \frac{\text{cm}^3}{\text{g}\cdot\text{s}}$.

**Parameter t<sub>d</sub>:** A study that infected human fibrosarcoma cells with Dengue Virus found an increase in phosphorylated eIF2α — indicating the activation of Integrated Stress Response arm of the UPR — at approximately 6 hours post infection [41]. At 24 – 36 hours in that same study, the inositol-requiring-protein-1 UPR pathway was activated. This suggests a timescale for the UPR may be in the loose range of 6 – 36 hours, and we use the mean of the extremes, $t_d \approx 21$ hours or $t_d \approx 80000$ s.

**Parameter A:** The most challenging parameter to estimate is the maximum flux of P-proteins. We first note that a heterodimer model on hamster data estimates the rate of monomer production to be $\lambda = 4 \frac{\mu\text{g}}{(\text{g brain})\cdot\text{day}}$ [40], which with brain density of 1 g/cm<sup>3</sup> [39], we have $\lambda \approx 5 \times 10^{-11} \frac{\text{g}}{\text{cm}^3\text{s}}$.

There are approximately $2 \times 10^9$ neurons per cm<sup>3</sup> [42] across many species. Neuronal perikaryal diameters range from about 6 to 80 µm and the perikarya accounts for approximately 10% of the neuron surface area [43]. By these estimate and with the geometric mean of the diameter lower and upper bounds taken as a length scale, the surface area of a neuron is $\sim 5000 \text{ µm}^2$ or $5 \times 10^{-5} \text{ cm}^2$.

Then, as a ballpark figure for the flux of P-proteins from a single neuron, we obtain $A^* \approx \frac{5 \times 10^{-11}}{2 \times 10^9 \times 5 \times 10^{-5}} = 5 \times 10^{-16} \frac{\text{g}}{\text{cm}^2 \cdot \text{s}}$. To this figure, we make a scaling correction. As our intent with the model is to understand the dynamics experienced by a single typical neuron, we need to take into account that no neuron functions in isolation and the brain is full of neurons. Neurons in a nearby vicinity will also experience similar stresses due to the presence of S-proteins. We thereby choose to scale $A^*$ by a ballpark estimate of the number of neurons that could influence any particular neuron through their release of P-proteins.

On average, a P-protein survives a time of $1/a$, resulting in a characteristic diffusive lengthscale of $\sqrt{D_p}/a \approx 0.06$ cm. Within a volume of $0.06^3 \text{ cm}^3$, there are approximately $4 \times 10^5$ neurons, and we scale $A^*$ up by this to give a nominal flux of $A \approx 2 \times 10^{-10} \frac{\text{g}}{\text{cm}^2 \cdot \text{s}}$.

**Parameter R:** cells of the globus pallidus in humans can have dendritic extends of $\approx 1000 \text{ µm}$ and soma radii of $\approx 16.5 \text{ µm}$ [44]. To come up with a representative length scale for an idealized “spherical neuron”, we take a geometric mean of these length scales giving $R \sim 128 \text{ µm} \approx 0.01$ cm.

**Setup for Systems** We focus on the behavior of our model with the estimated biological parameters. Numerical specifications for our biologically-driven study are given in Appendix B.2. The parameters of Table 4 define characteristic scales of $\bar{t} = 20000$ s, $\bar{x} = 0.06$ cm, $\bar{P} = 6 \times 10^{-5}$ g/cm<sup>3</sup>, and $\bar{S} = 0.0001$ g/cm<sup>3</sup>. Throughout all simulations, we choose initial conditions so that the
<table><caption>Table 4: Loose estimates of parameter values. *: value reported is 4 × 10⁵ times the true estimated flux of a single neuron in isolation — see explanation in text.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td>D<sub>P</sub></td><td>2 × 10<sup>-7</sup> cm<sup>2</sup>/s</td><td>[33, 34, 35]</td></tr><tr><td>D<sub>S</sub></td><td>2 × 10<sup>-8</sup> cm<sup>2</sup>/s</td><td>[33, 34, 37, 35]</td></tr><tr><td>a</td><td>5 × 10<sup>-5</sup> s<sup>-1</sup></td><td>[37]</td></tr><tr><td>b</td><td>9 × 10<sup>-7</sup> s<sup>-1</sup></td><td>[37]</td></tr><tr><td>S<sub>c</sub></td><td>1 × 10<sup>-4</sup> g/cm<sup>3</sup></td><td>[38, 39]</td></tr><tr><td>f</td><td>0 s<sup>-1</sup></td><td>N/A</td></tr><tr><td>c</td><td>2 cm<sup>3</sup>/s</td><td>[40, 39]</td></tr><tr><td>t<sub>d</sub></td><td>8 × 10<sup>4</sup> s</td><td>[41]</td></tr><tr><td>A</td><td>2 × 10<sup>-10</sup> g/cm<sup>2</sup>s</td><td>[40, 39, 42, 43]</td></tr><tr><td>R</td><td>1 × 10<sup>-2</sup> cm</td><td>[44]</td></tr></tbody></table>

<table><caption>Table 5: Dimensionless values for clinical model.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>γ</td><td>4</td></tr><tr><td>σ</td><td>0</td></tr><tr><td>δ</td><td>0.1</td></tr><tr><td>η</td><td>3.16</td></tr><tr><td>β</td><td>0.018</td></tr><tr><td>τ</td><td>4</td></tr><tr><td>ρ</td><td>0.158</td></tr></tbody></table>
<table><caption>Table 4: Loose estimates of parameter values. *: value reported is 4 &times; 10<sup>5</sup> times the true estimated flux of a single neuron in isolation — see explanation in text.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td>D<sub>P</sub></td><td>2 &times; 10<sup>-7</sup> cm<sup>2</sup>/s</td><td>[33, 34, 35]</td></tr><tr><td>D<sub>S</sub></td><td>2 &times; 10<sup>-8</sup> cm<sup>2</sup>/s</td><td>[33, 34, 37, 35]</td></tr><tr><td>a</td><td>5 &times; 10<sup>-5</sup> s<sup>-1</sup></td><td>[37]</td></tr><tr><td>b</td><td>9 &times; 10<sup>-7</sup> s<sup>-1</sup></td><td>[37]</td></tr><tr><td>S<sub>c</sub></td><td>1 &times; 10<sup>-4</sup> g/cm<sup>3</sup></td><td>[38, 39]</td></tr><tr><td>f</td><td>0 s<sup>-1</sup></td><td>N/A</td></tr><tr><td>c</td><td>2 cm<sup>3</sup>/s</td><td>[40, 39]</td></tr><tr><td>t<sub>d</sub></td><td>8 &times; 10<sup>4</sup> s</td><td>[41]</td></tr><tr><td>A</td><td>2 &times; 10<sup>-10</sup> g/cm<sup>2</sup>&middot;s</td><td>[40, 39, 42, 43]</td></tr><tr><td>R</td><td>1 &times; 10<sup>-2</sup> cm</td><td>[44]</td></tr></tbody></table>

<table><caption>Table 5: Dimensionless values for clinical model.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>&gamma;</td><td>4</td></tr><tr><td>&sigma;</td><td>0</td></tr><tr><td>&delta;</td><td>0.1</td></tr><tr><td>&eta;</td><td>3.16</td></tr><tr><td>&beta;</td><td>0.018</td></tr><tr><td>&tau;</td><td>4</td></tr><tr><td>&rho;</td><td>0.158</td></tr></tbody></table>
dimensional values of $i_c$ and $i_s$ of Eqs. (27) and (28) are

$$ I_P = \bar{P}\bar{x}^d \approx \begin{cases} 4 \times 10^{-6} \text{g/cm}^2, & d=1 \\ 3 \times 10^{-7} \text{g/cm}, & d=2 \\ 2 \times 10^{-8} \text{g}, & d=3 \end{cases} \qquad (34) $$

$$ I_S = \bar{S}\bar{x}^d \approx \begin{cases} 6 \times 10^{-6} \text{g/cm}^2, & d=1 \\ 4 \times 10^{-7} \text{g/cm}, & d=2 \\ 3 \times 10^{-8} \text{g}, & d=3 \end{cases} \qquad (35) $$

when $i_p = i_s = 1$. We also ensure the initial concentration profiles match exactly in dimensional units across all systems.

To retain a consistent dimensional point of comparison with identical initial conditions, we nondimensionalize according to the baselines parameters and simulate; this results in a slightly more general nondimensionalization discussed in Appendix A.1.

## 4.2 Sensitivity Analysis

The capacity for medicine to modify some or most of the parameters seems likely. We therefore consider how small, 30%, changes in each modifiable parameter in the system affect the UPR mechanism. Specifically, we look at $h_{\text{avg}}$ (Eq. (25)) and

$$ W_{\text{avg}} := \omega_{\text{avg}}/\bar{t}, \qquad (36) $$

the dimensional value of $\omega_{\text{avg}}$ (Eq. (26)). We comment that $h_{\text{avg}}$ is already dimensionless.

For baseline, we use the parameter values of Table 4. Then, we vary each parameter individually by 30% up and down, keeping all others fixed at baseline. We assume cell sizes are fixed and do not change $R$. The results are found in Tables 6 and 7 in differing dimensions.

In the case $d=1$, the model has large $h_{\text{avg}}$-values, suggesting a very diseased state. In that case, even 30% changes in many parameters do not result in large changes because so many parameters contribute to high S-presence. We notice much larger changes with $d=2$ where the S-presence is smaller. For $d=1$, the most beneficial changes were decreasing $A$, increasing $S_c$, increasing $b$, and increasing $D_S$. For $d=2$, the most beneficial changes were decreasing $A$, decreasing $c$, increasing $D_S$, and increasing $S_c$. Both $d=1$ and $d=2$ show the same pattern of increase/decrease with the different parameter changes. The $d=3$ cases did not show notable disease prevalence with $h_{\text{avg}}=0$ and no oscillations.

As for the oscillation frequency, $W_{\text{avg}}$, most of the parameters had little effect. The only parameter that had an appreciable effect was the UPR delay $t_d$. Oscillations tended not to occur much above $d=1$, but when present in $d=1$ and $d=2$, they tended to be very close in value.

## 4.3 Case Study

It appears the drug pentosan polysulfate (PPS) can extend survival times in patients with prion diseases [45]. One man survived 10 years with CJD while receiving that drug treatment [46]. Researchers have suggested the drug inhibiting the binding of $PrP^{Sc}$ to $PrP^C$ or causing the fragmentation of $PrP^{Sc}$ as possible mechanisms [47]. It appears PPS does not affect levels of $PrP^C$ [47].
dimensional values of $i_c$ and $i_s$ of Eqs. (27) and (28) are

$$ I_P = \bar{P}\bar{x}^d \approx \begin{cases} 4 \times 10^{-6}\text{g/cm}^2, & d=1 \\ 3 \times 10^{-7}\text{g/cm}, & d=2 \\ 2 \times 10^{-8}\text{g}, & d=3 \end{cases} \quad (34) $$

$$ I_S = \bar{S}\bar{x}^d \approx \begin{cases} 6 \times 10^{-6}\text{g/cm}^2, & d=1 \\ 4 \times 10^{-7}\text{g/cm}, & d=2 \\ 3 \times 10^{-8}\text{g}, & d=3 \end{cases} \quad (35) $$

when $i_p = i_s = 1$. We also ensure the initial concentration profiles match exactly in dimensional units across all systems.

To retain a consistent dimensional point of comparison with identical initial conditions, we nondimensionalize according to the baselines parameters and simulate; this results in a slightly more general nondimensionalization discussed in Appendix A.1.

## 4.2 Sensitivity Analysis

The capacity for medicine to modify some or most of the parameters seems likely. We therefore consider how small, 30%, changes in each modifiable parameter in the system affect the UPR mechanism. Specifically, we look at $h_{\text{avg}}$ (Eq. (25)) and

$$ W_{\text{avg}} := \omega_{\text{avg}}/\bar{t}, \quad (36) $$

the dimensional value of $\omega_{\text{avg}}$ (Eq. (26)). We comment that $h_{\text{avg}}$ is already dimensionless.

For baseline, we use the parameter values of Table 4. Then, we vary each parameter individually by 30% up and down, keeping all others fixed at baseline. We assume cell sizes are fixed and do not change $R$. The results are found in Tables 6 and 7 in differing dimensions.

In the case $d=1$, the model has large $h_{\text{avg}}$-values, suggesting a very diseased state. In that case, even 30% changes in many parameters do not result in large changes because so many parameters contribute to high S-presence. We notice much larger changes with $d=2$ where the S-presence is smaller. For $d=1$, the most beneficial changes were decreasing $A$, increasing $S_c$, increasing $b$, and increasing $D_S$. For $d=2$, the most beneficial changes were decreasing $A$, decreasing $c$, increasing $D_S$, and increasing $S_c$. Both $d=1$ and $d=2$ show the same pattern of increase/decrease with the different parameter changes. The $d=3$ cases did not show notable disease prevalence with $h_{\text{avg}}=0$ and no oscillations.

As for the oscillation frequency, $W_{\text{avg}}$, most of the parameters had little effect. The only parameter that had an appreciable effect was the UPR delay $t_d$. Oscillations tended not to occur much above $d=1$, but when present in $d=1$ and $d=2$, they tended to be very close in value.

## 4.3 Case Study

It appears the drug pentosan polysulfate (PPS) can extend survival times in patients with prion diseases [45]. One man survived 10 years with CJD while receiving that drug treatment [46]. Researchers have suggested the drug inhibiting the binding of $PrP^{Sc}$ to $PrP^C$ or causing the fragmentation of $PrP^{Sc}$ as possible mechanisms [47]. It appears PPS does not affect levels of $PrP^C$ [47].
<table><caption>Table 6: Values of $h_{\text{avg}}$ under various parameter changes with changes relative to baseline in parentheses. The $h_{\text{avg}}$ values were 0 for all 3D cases.</caption><thead><tr><th>Parameter Changes</th><th><i>h</i><sub>avg</sub> (1D)</th><th><i>h</i><sub>avg</sub> (2D)</th></tr></thead><tbody><tr><td>baseline</td><td>0.830</td><td>0.0687</td></tr><tr><td><i>D</i><sub>P</sub> ← 0.7<i>D</i><sub>P</sub></td><td><span style="color:red;">0.831</span> (+0.1%)</td><td><span style="color:red;">0.124</span> (+80.5%)</td></tr><tr><td><i>D</i><sub>P</sub> ← 1.3<i>D</i><sub>P</sub></td><td><span style="color:blue;">0.828</span> (-0.2%)</td><td><span style="color:blue;">0.032</span> (-53.7%)</td></tr><tr><td><i>D</i><sub>S</sub> ← 0.7<i>D</i><sub>S</sub></td><td><span style="color:red;">0.842</span> (+1.4%)</td><td><span style="color:red;">0.187</span> (+172%)</td></tr><tr><td><i>D</i><sub>S</sub> ← 1.3<i>D</i><sub>S</sub></td><td><span style="color:blue;">0.812</span> (-2.2%)</td><td><span style="color:blue;">0.0107</span> (-84.4%)</td></tr><tr><td><i>a</i> ← 0.7<i>a</i></td><td><span style="color:red;">0.840</span> (+1.2%)</td><td><span style="color:red;">0.108</span> (+57.2%)</td></tr><tr><td><i>a</i> ← 1.3<i>a</i></td><td><span style="color:blue;">0.821</span> (-1.1%)</td><td><span style="color:blue;">0.037</span> (-46.1%)</td></tr><tr><td><i>b</i> ← 0.7<b></i></td><td><span style="color:red;">0.847</span> (+2.0%)</td><td><span style="color:red;">0.0947</span> (+37.8%)</td></tr><tr><td><i>b</i> ← 1.3<b></i></td><td><span style="color:blue;">0.810</span> (-2.4%)</td><td><span style="color:blue;">0.047</span> (-31.6%)</td></tr><tr><td><i>c</i> ← 0.7<i>c</i></td><td><span style="color:blue;">0.815</span> (-1.8%)</td><td><span style="color:blue;">0.00211</span> (-96.9%)</td></tr><tr><td><i>c</i> ← 1.3<i>c</i></td><td><span style="color:red;">0.838</span> (+1.0%)</td><td><span style="color:red;">0.140</span> (+103.8%)</td></tr><tr><td><i>A</i> ← 0.7<i>A</i></td><td><span style="color:blue;">0.772</span> (-7.0%)</td><td><span style="color:blue;">0.000128</span> (-99.8%)</td></tr><tr><td><i>A</i> ← 1.3<i>A</i></td><td><span style="color:red;">0.860</span> (+3.6%)</td><td><span style="color:red;">0.218</span> (+217.3%)</td></tr><tr><td><i>t</i><sub>d</sub> ← 0.7<i>t</i><sub>d</sub></td><td><span style="color:red;">0.836</span> (+0.7%)</td><td><span style="color:red;">0.0696</span> (+1.3%)</td></tr><tr><td><i>t</i><sub>d</sub> ← 1.3<i>t</i><sub>d</sub></td><td><span style="color:blue;">0.823</span> (-0.8%)</td><td><span style="color:blue;">0.0678</span> (-13.5%)</td></tr><tr><td><i>S</i><sub>c</sub> ← 0.7<i>S</i><sub>c</sub></td><td><span style="color:red;">0.861</span> (+3.7%)</td><td><span style="color:red;">0.160</span> (+132.9%)</td></tr><tr><td><i>S</i><sub>c</sub> ← 1.3<i>S</i><sub>c</sub></td><td><span style="color:blue;">0.792</span> (-4.6%)</td><td><span style="color:blue;">0.0197</span> (-71.3%)</td></tr></tbody></table>

<table><caption>Table 7: Values of $W_{\text{avg}}$ at various parameter changes. No oscillations were noted for any 3D cases.</caption><thead><tr><th>Parameter Changes</th><th><i>W</i><sub>avg</sub> [μHz] (1D)</th><th><i>W</i><sub>avg</sub> [μHz] (2D)</th></tr></thead><tbody><tr><td>baseline</td><td>3.85</td><td>0</td></tr><tr><td><i>D</i><sub>P</sub> ← 0.7<i>D</i><sub>P</sub></td><td>3.91</td><td>0</td></tr><tr><td><i>D</i><sub>P</sub> ← 1.3<i>D</i><sub>P</sub></td><td>3.79</td><td>0</td></tr><tr><td><i>D</i><sub>S</sub> ← 0.7<i>D</i><sub>S</sub></td><td>3.79</td><td>4.06</td></tr><tr><td><i>D</i><sub>S</sub> ← 1.3<i>D</i><sub>S</sub></td><td>3.87</td><td>0</td></tr><tr><td><i>a</i> ← 0.7<i>a</i></td><td>3.84</td><td>0</td></tr><tr><td><i>a</i> ← 1.3<i>a</i></td><td>3.85</td><td>0</td></tr><tr><td><i>b</i> ← 0.7<b></i></td><td>3.78</td><td>0</td></tr><tr><td><i>b</i> ← 1.3<b></i></td><td>3.80</td><td>0</td></tr><tr><td><i>c</i> ← 0.7<b></i></td><td>3.74</td><td>0</td></tr><tr><td><i>c</i> ← 1.3<b></i></td><td>3.92</td><td>4.10</td></tr><tr><td><i>A</i> ← 0.7<i>A</i></td><td>3.81</td><td>0</td></tr><tr><td><i>A</i> ← 1.3<i>A</i></td><td>3.87</td><td>4.20</td></tr><tr><td><i>t</i><sub>d</sub> ← 0.7<i>t</i><sub>d</sub></td><td>5.19</td><td>0</td></tr><tr><td><i>t</i><sub>d</sub> ← 1.3<i>t</i><sub>d</sub></td><td>3.02</td><td>0</td></tr><tr><td><i>S</i><sub>c</sub> ← 0.7<i>S</i><sub>c</sub></td><td>3.76</td><td>4.03</td></tr><tr><td><i>S</i><sub>c</sub> ← 1.3<i>S</i><sub>c</sub></td><td>3.88</td><td>0</td></tr></tbody></table>
<table><caption>Table 6: Values of $h_{\text{avg}}$ under various parameter changes with changes relative to baseline in parentheses. The $h_{\text{avg}}$ values were 0 for all 3D cases.</caption><thead><tr><th>Parameter Changes</th><th><i>h</i><sub>avg</sub> (1D)</th><th><i>h</i><sub>avg</sub> (2D)</th></tr></thead><tbody><tr><td>baseline</td><td>0.830</td><td>0.0687</td></tr><tr><td><i>D</i><sub>P</sub> ← 0.7<i>D</i><sub>P</sub></td><td><span style="color:red;">0.831</span> (+0.1%)</td><td><span style="color:red;">0.124</span> (+80.5%)</td></tr><tr><td><i>D</i><sub>P</sub> ← 1.3<i>D</i><sub>P</sub></td><td><span style="color:blue;">0.828</span> (-0.2%)</td><td><span style="color:blue;">0.032</span> (-53.7%)</td></tr><tr><td><i>D</i><sub>S</sub> ← 0.7<i>D</i><sub>S</sub></td><td><span style="color:red;">0.842</span> (+1.4%)</td><td><span style="color:red;">0.187</span> (+172%)</td></tr><tr><td><i>D</i><sub>S</sub> ← 1.3<i>D</i><sub>S</sub></td><td><span style="color:blue;">0.812</span> (-2.2%)</td><td><span style="color:blue;">0.0107</span> (-84.4%)</td></tr><tr><td><i>a</i> ← 0.7<i>a</i></td><td><span style="color:red;">0.840</span> (+1.2%)</td><td><span style="color:red;">0.108</span> (+57.2%)</td></tr><tr><td><i>a</i> ← 1.3<i>a</i></td><td><span style="color:blue;">0.821</span> (-1.1%)</td><td><span style="color:blue;">0.037</span> (-46.1%)</td></tr><tr><td><i>b</i> ← 0.7<i>b</i></td><td><span style="color:red;">0.847</span> (+2.0%)</td><td><span style="color:red;">0.0947</span> (+37.8%)</td></tr><tr><td><i>b</i> ← 1.3<i>b</i></td><td><span style="color:blue;">0.810</span> (-2.4%)</td><td><span style="color:blue;">0.047</span> (-31.6%)</td></tr><tr><td><i>c</i> ← 0.7<i>c</i></td><td><span style="color:blue;">0.815</span> (-1.8%)</td><td><span style="color:blue;">0.00211</span> (-96.9%)</td></tr><tr><td><i>c</i> ← 1.3<i>c</i></td><td><span style="color:red;">0.838</span> (+1.0%)</td><td><span style="color:red;">0.140</span> (+103.8%)</td></tr><tr><td><i>A</i> ← 0.7<i>A</i></td><td><span style="color:blue;">0.772</span> (-7.0%)</td><td><span style="color:blue;">0.000128</span> (-99.8%)</td></tr><tr><td><i>A</i> ← 1.3<i>A</i></td><td><span style="color:red;">0.860</span> (+3.6%)</td><td><span style="color:red;">0.218</span> (+217.3%)</td></tr><tr><td><i>t</i><sub>d</sub> ← 0.7<i>t</i><sub>d</sub></td><td><span style="color:red;">0.836</span> (+0.7%)</td><td><span style="color:red;">0.0696</span> (+1.3%)</td></tr><tr><td><i>t</i><sub>d</sub> ← 1.3<i>t</i><sub>d</sub></td><td><span style="color:blue;">0.823</span> (-0.8%)</td><td><span style="color:blue;">0.0678</span> (-13.5%)</td></tr><tr><td><i>S</i><sub>c</sub> ← 0.7<i>S</i><sub>c</sub></td><td><span style="color:red;">0.861</span> (+3.7%)</td><td><span style="color:red;">0.160</span> (+132.9%)</td></tr><tr><td><i>S</i><sub>c</sub> ← 1.3<i>S</i><sub>c</sub></td><td><span style="color:blue;">0.792</span> (-4.6%)</td><td><span style="color:blue;">0.0197</span> (-71.3%)</td></tr></tbody></table>

<table><caption>Table 7: Values of $W_{\text{avg}}$ at various parameter changes. No oscillations were noted for any 3D cases.</caption><thead><tr><th>Parameter Changes</th><th><i>W</i><sub>avg</sub> [μHz] (1D)</th><th><i>W</i><sub>avg</sub> [μHz] (2D)</th></tr></thead><tbody><tr><td>baseline</td><td>3.85</td><td>0</td></tr><tr><td><i>D</i><sub>P</sub> ← 0.7<i>D</i><sub>P</sub></td><td>3.91</td><td>0</td></tr><tr><td><i>D</i><sub>P</sub> ← 1.3<i>D</i><sub>P</sub></td><td>3.79</td><td>0</td></tr><tr><td><i>D</i><sub>S</sub> ← 0.7<i>D</i><sub>S</sub></td><td>3.79</td><td>4.06</td></tr><tr><td><i>D</i><sub>S</sub> ← 1.3<i>D</i><sub>S</sub></td><td>3.87</td><td>0</td></tr><tr><td><i>a</i> ← 0.7<i>a</i></td><td>3.84</td><td>0</td></tr><tr><td><i>a</i> ← 1.3<i>a</i></td><td>3.85</td><td>0</td></tr><tr><td><i>b</i> ← 0.7<i>b</i></td><td>3.78</td><td>0</td></tr><tr><td><i>b</i> ← 1.3<i>b</i></td><td>3.80</td><td>0</td></tr><tr><td><i>c</i> ← 0.7<i>c</i></td><td>3.74</td><td>0</td></tr><tr><td><i>c</i> ← 1.3<i>c</i></td><td>3.92</td><td>4.10</td></tr><tr><td><i>A</i> ← 0.7<i>A</i></td><td>3.81</td><td>0</td></tr><tr><td><i>A</i> ← 1.3<i>A</i></td><td>3.87</td><td>4.20</td></tr><tr><td><i>t</i><sub>d</sub> ← 0.7<i>t</i><sub>d</sub></td><td>5.19</td><td>0</td></tr><tr><td><i>t</i><sub>d</sub> ← 1.3<i>t</i><sub>d</sub></td><td>3.02</td><td>0</td></tr><tr><td><i>S</i><sub>c</sub> ← 0.7<i>S</i><sub>c</sub></td><td>3.76</td><td>4.03</td></tr><tr><td><i>S</i><sub>c</sub> ← 1.3<i>S</i><sub>c</sub></td><td>3.88</td><td>0</td></tr></tbody></table>
<figure><img src="image_47.png" /><figcaption>Figure 12: Variations in $h_{\text{avg}}$ and $W_{\text{avg}}$ at various potencies.</figcaption></figure>

Other drugs have also been studied where $PrP^C$ levels were made significantly lower — through chloropromazine (CPZ) and U18666A [47].

To speculate quantitatively on drug treatments in our model, we consider a hypothetical concoction of drugs that (1) increases the clearance rate of $PrP^{Sc}$, (2) increases the diffusivity of $PrP^{Sc}$, (3) reduces the capacity of $PrP^{Sc}$ to convert $PrP^C$, and (4) reduces the production of $PrP^C$. Items (1)-(2) are inspired by the speculation on the PPS drug, where we assume fragmented $PrP^{Sc}$ will be cleared more readily and diffuse more easily. Item (3) is again a possible benefit of PPS, and item (4) is based on the possible effects of CPZ and U18666A. Numerically, we investigate how $h_{\text{avg}}$ and $W_{\text{avg}}$ change as the drug “potency” changes. We define the potency $\lambda \in [1, \infty)$ so that at potency $\lambda$, the values of $D_s$ and $b$ are both increased by a factor $\lambda$ relative to baseline (Table 4) and the values of $c$ and $A$ are both decreased by a factor of $\lambda$ relative to baseline (Table 4). The resulting outcomes are found in Figure 12. We observe that $h_{\text{avg}}$ is nearly zero and oscillations stop above $\lambda \approx 2.4$.

## 4.4 Clinical Relevance

At this point, we discuss the potential clinical significance of our model. From Table 6, we note that the change in each dimensional parameter has a corresponding change in the strength of the UPR response $h_{\text{avg}}$. This suggests that medical interventions that increase the diffusivities of either or both proteins; reduce the P-S conversion rate; enhance the rate that either or both proteins are cleared; shorten the UPR delay; decrease the maximum flux of P; and increase the threshold at which the UPR is triggered, may benefit patients clinically. The most significant reductions in the UPR intensity seem to come about through decreasing the maximum rate P is released, decreasing the recruitment rate, increasing the diffusivity of S, increasing the rate at which S is cleared, and increasing the threshold sensitivity of the UPR to S. We note these results may vary depending on the point about which sensitivity is studied. We remark that based on the possible effects of drugs used in prion disease treatment to date, their methods of action are strongly aligned with these targets.

An experimentalist may be able to identify signs of the UPR mechanism through reduced production of $PrP^C$. Based on our simulations, we expect that oscillations will become more pronounced when the relative diffusivity of $PrP^{Sc}$ relative to $PrP^C$ is small, the relative clearance rate of $PrP^{Sc}$ relative to $PrP^C$ is small, and the UPR delay is large relative to the $PrP^C$ clearance time. It is also expected that the temporal frequency of oscillations will decrease as the UPR delay time increases
<figure><img src="image_45.png" /><figcaption>Figure 12: Variations in $h_{\text{avg}}$ and $W_{\text{avg}}$ at various potencies.</figcaption></figure>

Other drugs have also been studied where $PrP^C$ levels were made significantly lower — through chloropromazine (CPZ) and U18666A [47].

To speculate quantitatively on drug treatments in our model, we consider a hypothetical con-
cption of drugs that (1) increases the clearance rate of $PrP^{Sc}$, (2) increases the diffusivity of $PrP^{Sc}$,
(3) reduces the capacity of $PrP^{Sc}$ to convert $PrP^C$, and (4) reduces the production of $PrP^C$. Items
(1)-(2) are inspired by the speculation on the PPS drug, where we assume fragmented $PrP^{Sc}$ will be
cleared more readily and diffuse more easily. Item (3) is again a possible benefit of PPS, and item
(4) is based on the possible effects of CPZ and U18666A. Numerically, we investigate how $h_{\text{avg}}$ and
$W_{\text{avg}}$ change as the drug “potency” changes. We define the potency $\lambda \in [1, \infty)$ so that at potency
$\lambda$, the values of $D_s$ and $b$ are both increased by a factor $\lambda$ relative to baseline (Table 4) and the
values of $c$ and $A$ are both decreased by a factor of $\lambda$ relative to baseline (Table 4). The resulting
outcomes are found in Figure 12. We observe that $h_{\text{avg}}$ is nearly zero and oscillations stop above
$\lambda \approx 2.4$.

## 4.4 Clinical Relevance

At this point, we discuss the potential clinical significance of our model. From Table 6, we note
that the change in each dimensional parameter has a corresponding change in the strength of the
UPR response $h_{\text{avg}}$. This suggests that medical interventions that increase the diffusivities of either
or both proteins; reduce the P-S conversion rate; enhance the rate that either or both proteins are
cleared; shorten the UPR delay; decrease the maximum flux of P; and increase the threshold at
which the UPR is triggered, may benefit patients clinically. The most significant reductions in the
UPR intensity seem to come about through decreasing the maximum rate P is released, decreasing
the recruitment rate, increasing the diffusivity of S, increasing the rate at which S is cleared, and
increasing the threshold sensitivity of the UPR to S. We note these results may vary depending on
the point about which sensitivity is studied. We remark that based on the possible effects of drugs
used in prion disease treatment to date, their methods of action are strongly aligned with these
targets.

An experimentalist may be able to identify signs of the UPR mechanism through reduced produc-
tion of $PrP^C$. Based on our simulations, we expect that oscillations will become more pronounced
when the relative diffusivity of $PrP^{Sc}$ relative to $PrP^C$ is small, the relative clearance rate of $PrP^{Sc}$
relative to $PrP^C$ is small, and the UPR delay is large relative to the $PrP^C$ clearance time. It is also
expected that the temporal frequency of oscillations will decrease as the UPR delay time increases
and that for large enough diffusivities or clearances of $PrP^{Sc}$, the oscillations may entirely fade away.

## 4.5 Model Limitations

To arrive at the system we studied, we had to make several assumptions. A proper understanding of the resulting limitations is indispensable for healthy scientific scepticism in interpreting the model results and in developing better, more insightful models. We include a list of limitations below along with possible remedies.

* **Only two proteins:** instead of including multitudes of oligomeric states and stable nuclei, we chose to only include P- and S-proteins. A more refined model could consider assortments of assembly sizes, with associated mass-dependent diffusivities and clearance rates.

* **Only one neuron:** we centered our model on one neuron, but it is possible to consider multiple neurons, with associated positions in space so that species that diffuse away from one neuron can affect another.

* **Simplified neuron geometry:** we used simplified geometries for the neuron shape. Due to the long axons relative to the cell body size, radial symmetry is not accurate. Through a finite element approach, a more complicated geometry could be tackled.

* **Constant reaction rates:** through the ageing process, the rates of misfolding or clearance, among other parameters, may vary. It is possible to consider time-dependent parameters as well.

* **No membrane-bound P-protein:** we did not consider a form of P-protein that is bound to the cell membrane. Further refinements could add this.

* **Protein production not impacted by UPR:** the model could be adapted to include a reduction in maximal P-protein production due to cumulative effects of the UPR. Likewise, since the UPR may be a mechanism of cell death, the model could be extended to determine when the neuron dies.

* **Only one cell type:** there are many cell types in the brain, but we focused only on neurons. By including other cell types, we could better model the effects of stress and inflammation.

* **Lack of data-driven parameters:** due to the uncertainty of the different model parameters, we focused mostly on the qualitative aspects of the model. Our clinical exploration study provided interesting insights with regards to possible medical treatments, but accurate measurements of parameters are needed.

However, even with these limitations, our simple approach already highlights some important features of these complex mechanisms, as well as the role of key parameters involved.

## 5 Conclusion and Future Work

In this paper, we presented a nonlinear, coupled system of reaction-diffusion equations with non-linear, delayed boundary conditions to model the Unfolded Protein Response in a simplified setting
and that for large enough diffusivities or clearances of $PrP^{Sc}$, the oscillations may entirely fade away.

## 4.5 Model Limitations

To arrive at the system we studied, we had to make several assumptions. A proper understanding of the resulting limitations is indispensable for healthy scientific scepticism in interpreting the model results and in developing better, more insightful models. We include a list of limitations below along with possible remedies.

* **Only two proteins:** instead of including multitudes of oligomeric states and stable nuclei, we chose to only include P- and S-proteins. A more refined model could consider assortments of assembly sizes, with associated mass-dependent diffusivities and clearance rates.

* **Only one neuron:** we centered our model on one neuron, but it is possible to consider multiple neurons, with associated positions in space so that species that diffuse away from one neuron can affect another.

* **Simplified neuron geometry:** we used simplified geometries for the neuron shape. Due to the long axons relative to the cell body size, radial symmetry is not accurate. Through a finite element approach, a more complicated geometry could be tackled.

* **Constant reaction rates:** through the ageing process, the rates of misfolding or clearance, among other parameters, may vary. It is possible to consider time-dependent parameters as well.

* **No membrane-bound P-protein:** we did not consider a form of P-protein that is bound to the cell membrane. Further refinements could add this.

* **Protein production not impacted by UPR:** the model could be adapted to include a reduction in maximal P-protein production due to cumulative effects of the UPR. Likewise, since the UPR may be a mechanism of cell death, the model could be extended to determine when the neuron dies.

* **Only one cell type:** there are many cell types in the brain, but we focused only on neurons. By including other cell types, we could better model the effects of stress and inflammation.

* **Lack of data-driven parameters:** due to the uncertainty of the different model parameters, we focused mostly on the qualitative aspects of the model. Our clinical exploration study provided interesting insights with regards to possible medical treatments, but accurate measurements of parameters are needed.

However, even with these limitations, our simple approach already highlights some important features of these complex mechanisms, as well as the role of key parameters involved.

## 5 Conclusion and Future Work

In this paper, we presented a nonlinear, coupled system of reaction-diffusion equations with non-linear, delayed boundary conditions to model the Unfolded Protein Response in a simplified setting
with representative healthy and unhealthy proteins. We found that oscillations in neuronal activity may be found under certain parameters and that through modifying certain biological parameters it may be possible to lessen intensity of the UPR.

To extend our work and make it more applicable to neurological systems and the study of neurodegenerative diseases, we consider obtaining accurate estimates for the model parameters based on clinical and experimental data, incorporating additional biological features and realism (see limitations section), and coming to a theoretical understanding of the effects of each model parameter.

# Acknowledgments

This project has received support from NSF Award DMS #2150478 and NSF Award DMS #2316952. This project has received support from Agence National de la Recherche PrionDiff 350 Project-ANR-21-CE15-0011. The authors would like to thank Erwin Suazo and Tamer Oraby for funding students in this work.

# A Numerical Method

## A.1 Different Nondimensionalization

With scales $\bar{t}$, $\bar{x}$, $\bar{P}$, and $\bar{S}$ chosen from the baseline parameters, the dimensionless system, more generally, takes the form

$$p,t = \delta_p(p,rr + \frac{d-1}{r}p,r) - \gamma ps - (\alpha + \sigma)p, \quad (37)$$

$$s,t = \delta_s(s,rr + \frac{d-1}{r}s,r) + \eta\gamma ps + \eta\sigma p - \beta s, \quad (38)$$

$$p,r|_{r=\rho} = \frac{-\theta}{1 + (\mu s(\rho, t - \tau))^m} \quad (39)$$

$$h_{\text{avg}} = \lim_{\tau \to \infty} \int_{0}^{\tau} \left(1 - \frac{1}{1 + (\mu s(\rho, t - \tau))^m}\right) dt, \quad (40)$$

where we have suppressed equations that do not change their form. Here, $\delta_p = \frac{\bar{D}_P}{\bar{x}^2}$, $\gamma = \bar{S}\bar{t}c$, $\alpha = \bar{t}a$, $\sigma = \bar{t}f$, $\eta = \frac{\bar{P}}{\bar{S}}$, $\delta_s = \frac{\bar{D}_S}{\bar{x}^2}$, $\beta = \bar{t}b$, $\rho = \frac{\bar{R}}{\bar{x}}$, $\theta = \frac{\bar{x}A}{\bar{P}D_P}$, $\mu = \frac{\bar{S}}{\bar{S}_c}$, and $\tau = \frac{t_d}{t}$. With the baseline parameters, $\delta_p = \alpha = \theta = \mu = 1$.

## A.2 Implementation

We solve the dimensionless system of Appendix A.1. The domain for $r$, $[\rho, \infty)$, is unbounded and we need to truncate the computational domain at some value, $r_\infty$. From a back-of-the-envelope calculation, we note that for $s = 0$, a one-dimensional steady state solution is $p(r) = \frac{1}{\sqrt{1+\sigma}} \exp(-\sqrt{1+\sigma}(r-\rho)) \le \exp(-(r-\rho))$. And if we impose $p(r_\infty) = 0$ with $r_\infty = \rho+7$, the truncation error is $\sim 0.0009$. And $p$ and $s$ are coupled together, with $p$ being the only source for $s$. Solutions tend to decay faster in 2 and 3 dimensions. Thus, for $O(1)$ solutions and data, choosing $r_\infty = \rho+7$, we anticipate associated truncation errors to be $\sim 0.001$ in the far field.
with representative healthy and unhealthy proteins. We found that oscillations in neuronal activity may be found under certain parameters and that through modifying certain biological parameters it may be possible to lessen intensity of the UPR.

To extend our work and make it more applicable to neurological systems and the study of neurodegenerative diseases, we consider obtaining accurate estimates for the model parameters based on clinical and experimental data, incorporating additional biological features and realism (see limitations section), and coming to a theoretical understanding of the effects of each model parameter.

# Acknowledgments

This project has received support from NSF Award DMS #2150478 and NSF Award DMS #2316952. This project has received support from Agence National de la Recherche PrionDiff 350 Project-ANR-21-CE15-0011. The authors would like to thank Erwin Suazo and Tamer Oraby for funding students in this work.

# A Numerical Method

## A.1 Different Nondimensionalization

With scales $\bar{t}$, $\bar{x}$, $\bar{P}$, and $\bar{S}$ chosen from the baseline parameters, the dimensionless system, more generally, takes the form

$$p,t = \delta_p(p,rr + \frac{d-1}{r}p,r) - \gamma ps - (\alpha + \sigma)p, \quad (37)$$

$$s,t = \delta_s(s,rr + \frac{d-1}{r}s,r) + \eta\gamma ps + \eta\sigma p - \beta s, \quad (38)$$

$$p,r|_{r=\rho} = \frac{-\theta}{1 + (\mu s(\rho, t - \tau))^m} \quad (39)$$

$$h_{\text{avg}} = \lim_{\tau \to \infty} \int_{0}^{\tau} \left(1 - \frac{1}{1 + (\mu s(\rho, t - \tau))^m}\right) dt, \quad (40)$$

where we have suppressed equations that do not change their form. Here, $\delta_p = \frac{\bar{D}_P}{\bar{x}^2}$, $\gamma = \bar{S}\bar{t}c$, $\alpha = \bar{t}a$, $\sigma = \bar{t}f$, $\eta = \frac{\bar{P}}{\bar{S}}$, $\delta_s = \frac{\bar{D}_S}{\bar{x}^2}$, $\beta = \bar{t}b$, $\rho = \frac{\bar{R}}{\bar{x}}$, $\theta = \frac{\bar{x}A}{\bar{P}D_P}$, $\mu = \frac{\bar{S}}{\bar{S}_c}$, and $\tau = \frac{t_d}{t}$. With the baseline parameters, $\delta_p = \alpha = \theta = \mu = 1$.

## A.2 Implementation

We solve the dimensionless system of Appendix A.1. The domain for $r$, $[\rho, \infty)$, is unbounded and we need to truncate the computational domain at some value, $r_\infty$. From a back-of-the-envelope calculation, we note that for $s = 0$, a one-dimensional steady state solution is $p(r) = \frac{1}{\sqrt{1+\sigma}} \exp(-\sqrt{1+\sigma}(r-\rho)) \le \exp(-(r-\rho))$. And if we impose $p(r_\infty) = 0$ with $r_\infty = \rho + 7$, the truncation error is $\sim 0.0009$. And $p$ and $s$ are coupled together, with $p$ being the only source for $s$. Solutions tend to decay faster in 2 and 3 dimensions. Thus, for $O(1)$ solutions and data, choosing $r_\infty = \rho + 7$, we anticipate associated truncation errors to be $\sim 0.001$ in the far field.
We choose positive integers $N_r$ and $N_t$ to be the meshing parameters in space and time. For a computational domain running from $r = \rho$ to $r = r_\infty$ and $t = 0$ to $t = t_\infty$, we define the spatial and temporal step sizes by

$$ \Delta r = \frac{r_\infty - \rho}{N_r}, \qquad \Delta t = \frac{t_\infty}{N_t}. $$

Then for $i = 0, 1, 2, \dots, N_r$ and $j = 0, 1, \dots, N_t$, we define $r_i = \rho + i\Delta r$ and $t_j = j\Delta t$. We denote $p_i^j$ as the numerical approximation to $p(x_i, t_j)$ and $s_i^j$ as the numerical approximation to $s(r_i, t_j)$. We furthermore take the delay $\tau = k\Delta t$ for some $k \in \{0, 1, 2, \dots\}$. At step $j$, second-order spatial discretizations of the diffusion and reaction terms for $p$ and $s$ are given by

$$ \bar{M}_p = \delta_p \frac{p_{i+1}^{j+1} - 2p_i^{j+1} + p_{i-1}^{j+1}}{\Delta r^2} - \gamma p_i^{j+1} s_i^j \\ - (1 + \sigma) p_i^{j+1}, \quad 0 \le i \le N_r - 1, \quad j = 0, 1, 2, \dots, \qquad (41) $$

$$ \bar{M}_s = \delta_s \frac{s_{i+1}^{j+1} - 2s_i^{j+1} + s_{i-1}^{j+1}}{\Delta r^2} + \eta\gamma p_i^j s_i^{j+1} \\ + \eta\sigma p_i^j - \beta s_i^{j+1}, \quad 0 \le i \le N_r - 1, \quad j = 0, 1, 2, \dots, \qquad (42) $$

respectively.

Handling the derivative conditions at $i = 0$ requires ghost points at the fictitious position $r_{-1} = \rho - \Delta r$, and the far-field at $i = N_r$, zero values are imposed. Thus:

$$ p_{N_r}^j = s_{N_r}^j = 0, \quad j = 0, 1, 2, \dots \qquad (43) $$

$$ \frac{p_1^j - p_{-1}^j}{2\Delta r} = \frac{-\theta}{1 + (\mu s_0^{j-k})^m}, \quad j = 0, 1, 2, \dots \qquad (44) $$

$$ \frac{s_1^j - s_{-1}^j}{2\Delta r} = 0, \quad j = 0, 1, 2, \dots \qquad (45) $$

To handle the delay, we assume $s_0^{-1}, s_0^{-2}, \dots, s_0^{-k}$ are given where $\tau = k\Delta t$. For initial conditions, we assume $p_i^0, s_i^0$ are given for $0 \le i \le N_r$.

These equations can be converted to matrix-vector form. We denote $p^j = (p_0^j, p_1^j, \dots, p_{N_r-1}^j)^T$ and $s^j = (s_0^j, s_1^j, \dots, s_{N_r-1}^j)^T$. Let $I \in \mathbb{R}^{N_r \times N_r}$ be the identity matrix, $D_2 \in \mathbb{R}^{N_r \times N_r}$ be given by

$$ D_2 = \frac{1}{\Delta r^2} \begin{pmatrix} -2 & 2 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & \vdots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 & -2 & 1 \\ 0 & 0 & \cdots & 0 & 1 & -2 \end{pmatrix}, $$
We choose positive integers $N_r$ and $N_t$ to be the meshing parameters in space and time. For a computational domain running from $r = \rho$ to $r = r_\infty$ and $t = 0$ to $t = t_\infty$, we define the spatial and temporal step sizes by

$$ \Delta r = \frac{r_\infty - \rho}{N_r}, \qquad \Delta t = \frac{t_\infty}{N_t}. $$

Then for $i = 0, 1, 2, \dots, N_r$ and $j = 0, 1, \dots, N_t$, we define $r_i = \rho + i\Delta r$ and $t_j = j\Delta t$. We denote $p_i^j$ as the numerical approximation to $p(x_i, t_j)$ and $s_i^j$ as the numerical approximation to $s(r_i, t_j)$. We furthermore take the delay $\tau = k\Delta t$ for some $k \in \{0, 1, 2, \dots\}$. At step $j$, second-order spatial discretizations of the diffusion and reaction terms for $p$ and $s$ are given by

$$ \bar{M}_p = \delta_p \frac{p_{i+1}^{j+1} - 2p_i^{j+1} + p_{i-1}^{j+1}}{\Delta r^2} - \gamma p_i^{j+1} s_i^j - (1 + \sigma) p_i^{j+1}, \quad 0 \le i \le N_r - 1, \quad j = 0, 1, 2, \dots, \qquad (41) $$

$$ \bar{M}_s = \delta_s \frac{s_{i+1}^{j+1} - 2s_i^{j+1} + s_{i-1}^{j+1}}{\Delta r^2} + \eta \gamma p_i^j s_i^{j+1} + \eta \sigma p_i^j - \beta s_i^{j+1}, \quad 0 \le i \le N_r - 1, \quad j = 0, 1, 2, \dots, \qquad (42) $$

respectively.

Handling the derivative conditions at $i = 0$ requires ghost points at the fictitious position $r_{-1} = \rho - \Delta r$, and the far-field at $i = N_r$, zero values are imposed. Thus:

$$ p_{N_r}^j = s_{N_r}^j = 0, \quad j = 0, 1, 2, \dots \qquad (43) $$

$$ \frac{p_1^j - p_{-1}^j}{2\Delta r} = \frac{-\theta}{1 + (\mu s_0^{j-k})^m}, \quad j = 0, 1, 2, \dots \qquad (44) $$

$$ \frac{s_1^j - s_{-1}^j}{2\Delta r} = 0, \quad j = 0, 1, 2, \dots \qquad (45) $$

To handle the delay, we assume $s_0^{-1}, s_0^{-2}, \dots, s_0^{-k}$ are given where $\tau = k\Delta t$. For initial conditions, we assume $p_i^0, s_i^0$ are given for $0 \le i \le N_r$.

These equations can be converted to matrix-vector form. We denote $p^j = (p_0^j, p_1^j, \dots, p_{N_r-1}^j)^\mathrm{T}$ and $s^j = (s_0^j, s_1^j, \dots, s_{N_r-1}^j)^\mathrm{T}$. Let $I \in \mathbb{R}^{N_r \times N_r}$ be the identity matrix, $D_2 \in \mathbb{R}^{N_r \times N_r}$ be given by

$$ D_2 = \frac{1}{\Delta r^2} \begin{pmatrix} -2 & 2 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & \vdots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 & -2 & 1 \\ 0 & 0 & \cdots & 0 & 1 & -2 \end{pmatrix}, $$
$D_1 \in \mathbb{R}^{N_r \times N_r}$ be given by

$$D_1 = \begin{cases} \begin{pmatrix} \frac{(d-1)}{2\Delta r} & 0 & 0 & \cdots & 0 & 0 \\ 0 & -1/r_1 & 0 & 1/r_1 & \cdots & 0 \\ 0 & \cdots & \cdots & \cdots & \cdots & \vdots \\ \vdots & \cdots & \cdots & \cdots & \ddots & \vdots \\ 0 & 0 & \cdots & -1/r_{N_r-2} & 0 & 1/r_{N_r-2} \\ 0 & 0 & \cdots & 0 & -1/r_{N_r-1} & 0 \end{pmatrix}, & d=2,3, \\ 0, & d=1, \end{cases}$$

and $V^j \in \mathbb{R}^{N_r \times 1}$ be given by

$$V^j = \left( \frac{2\theta\delta_p}{\Delta r(1 + (\mu s_0^{j-k})^m)}, 0, 0, \dots, 0 \right)^\top + \begin{cases} \left( \frac{-\theta\delta_p(d-1)}{\rho(1+(\mu s_0^{j-k})^m)}, 0, 0, \dots, 0 \right)^\top, & d=2,3, \\ 0, & d=1. \end{cases}$$

Defining functions $M_p, M_s : \mathbb{R}^{N_r \times 1} \to \mathbb{R}^{N_r \times N_r}$ via

$$M_p(s) = \delta_p(D_2 + D_1) - \gamma \operatorname{diag}(s) - (\sigma + \alpha)I \\ M_s(p) = \delta_s(D_2 + D_1) + \eta\gamma \operatorname{diag}(p) - \beta I,$$

we can define a weighted semi-implicit first-order time-stepping via

$$\begin{aligned} \frac{p^{j+1} - p^j}{\Delta t} &= w(M_p(s^j)p^{j+1} + V^{j+1}) + (1-w)(M_p(s^j)p^j + V^j) \\ \frac{s^{j+1} - s^j}{\Delta t} &= wM_s(p^j)s^{j+1} + (1-w)M_s(p^j)s^j + \eta\sigma p^j \end{aligned}$$

where $0 \le w \le 1$. With $w = 0$ this is the Euler method. With $w = 1$, this is a standard semi-implicit method. We take inspiration from the implicit Crank-Nicholson method [48] and take $w = 1/2$. In practice, we find $w = 1/2$ yields smaller errors with our semi-implicit approach. The semi-implicit systems require simple linear solves of two systems of size $N_r$. In contrast, Crank-Nicholson would require nonlinear solves of systems of size $2N_r$ and we find that the semi-implicit method is faster. When $k = 0$ (no delay), we compute $V^{j+1}$ by first estimating $s^{j+1}$ with the Improved Euler method.

Overall our method, as subsequently validated in Appendix A.4, is second order in space and first order in time.

## A.3 Frequency

To estimate the mean frequency over a finite computational window, we make some modifications. Because $s$ can take a while to peak and a final peak may take place sufficiently before $t_\infty$ to contribute a non-negligible error, we seek to estimate the mean peak-to-peak frequency of $s(\rho, t-\tau)$. But because oscillations may stop altogether, we need to avoid giving a frequency for situations where $s(\rho, t-\tau)$ peaks a few times early on and then stops. If there $N \le 1$ peaks for $s(\rho, t-\tau)$ on $(0, t_\infty)$, we define $\omega_{\text{avg}} = 0$ since a mean peak-to-peak interval cannot be computed for the frequency.
$D_1 \in \mathbb{R}^{N_r \times N_r}$ be given by

$$ D_1 = \left\{ \begin{array}{c} \frac{(d-1)}{2\Delta r} \\ \begin{array}{ccccccc} 0 & 0 & 0 & \cdots & 0 & 0 \\ -1/r_1 & 0 & 1/r_1 & \cdots & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & \vdots & \vdots \\ \vdots & \cdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & -1/r_{N_r-2} & 0 & 1/r_{N_r-2} \\ 0 & 0 & \cdots & 0 & -1/r_{N_r-1} & 0 \end{array} \\ 0, \quad d=1, \end{array} \right\} , \quad d=2,3, $$

and $V^j \in \mathbb{R}^{N_r \times 1}$ be given by

$$ V^j = \left( \frac{2\theta\delta_p}{\Delta r(1 + (\mu s_0^{j-k})^m)}, 0, 0, \dots, 0 \right)^\top + \begin{cases} \left( \frac{-\theta\delta_p(d-1)}{\rho(1+(\mu s_0^{j-k})^m)}, 0, 0, \dots, 0 \right)^\top, & d=2,3, \\ 0, & d=1. \end{cases} $$

Defining functions $M_p, M_s : \mathbb{R}^{N_r \times 1} \to \mathbb{R}^{N_r \times N_r}$ via

$$ M_p(s) = \delta_p(D_2 + D_1) - \gamma \operatorname{diag}(s) - (\sigma + \alpha)I $$

$$ M_s(p) = \delta_s(D_2 + D_1) + \eta\gamma \operatorname{diag}(p) - \beta I, $$

we can define a weighted semi-implicit first-order time-stepping via

$$ \frac{p^{j+1} - p^j}{\Delta t} = w(M_p(s^j)p^{j+1} + V^{j+1}) + (1-w)(M_p(s^j)p^j + V^j) $$

$$ \frac{s^{j+1} - s^j}{\Delta t} = w M_s(p^j) s^{j+1} + (1-w) M_s(p^j) s^j + \eta \sigma p^j $$

where $0 \le w \le 1$. With $w = 0$ this is the Euler method. With $w = 1$, this is a standard semi-implicit method. We take inspiration from the implicit Crank-Nicholson method [48] and take $w = 1/2$. In practice, we find $w = 1/2$ yields smaller errors with our semi-implicit approach. The semi-implicit systems require simple linear solves of two systems of size $N_r$. In contrast, Crank-Nicholson would require nonlinear solves of systems of size $2N_r$ and we find that the semi-implicit method is faster. When $k = 0$ (no delay), we compute $V^{j+1}$ by first estimating $s^{j+1}$ with the Improved Euler method.

Overall our method, as subsequently validated in Appendix A.4, is second order in space and first order in time.

## A.3 Frequency

To estimate the mean frequency over a finite computational window, we make some modifications. Because $s$ can take a while to peak and a final peak may take place sufficiently before $t_\infty$ to contribute a non-negligible error, we seek to estimate the mean peak-to-peak frequency of $s(\rho, t-\tau)$. But because oscillations may stop altogether, we need to avoid giving a frequency for situations where $s(\rho, t-\tau)$ peaks a few times early on and then stops. If there $N \le 1$ peaks for $s(\rho, t-\tau)$ on $(0, t_\infty)$, we define $\omega_{\text{avg}} = 0$ since a mean peak-to-peak interval cannot be computed for the frequency.
Assuming there are $N \ge 2$ peaks, let them occur at the times $0 < t_1 < t_2 < \dots < t_N < t_\infty$. We define a proxy for the mean frequency as

$$ \hat{\omega} = \frac{N-1}{t_N - t_1}. $$

Then we define

$$ \omega_{\text{avg}} = \begin{cases} \hat{\omega}, & t_{\infty} \le t_N + \max_{1 \le i \le N-1} (t_{i+1} - t_i) \\ 0, & \text{otherwise.} \end{cases} \qquad (46) $$

Intuitively, if the system is continuing to oscillate, then the end of the computational time $t_\infty$, should occur before the estimated latest possible time for the $(N+1)$st peak (otherwise there should have been $N+1$ peaks). If this condition fails, we define the frequency to be 0, indicating the system stopped oscillating.

## A.4 Validation of Numerical Scheme

Throughout the validation and simulations, the initial conditions prescribed are those of Equations (30)-(32).

To validate the convergence rates, we shall use mesh refinements. Let $N_r$ and $N_t$ be the meshing parameters. Denote $Q(N_r, N_t)_i^j = (p_i^j, q_i^j)^\top$ for $0 \le i \le N_r - 1, 0 \le j \le N_t$ to be the numerical approximation at $r = r_i$ and $t = t_j$ over the chosen grid. We define the operations

$$ \Lambda_r Q(N_r, N_t) = \max_{\substack{0 \le i \le N_r - 1 \\ 0 \le j \le N_t}} \|Q(N_r, N_t)_i^j - Q(2N_r, N_t)_{2i}^j\|_{\infty} $$

$$ \Lambda_t Q(N_r, N_t) = \max_{\substack{0 \le i \le N_r - 1 \\ 0 \le j \le N_t}} \|Q(N_r, N_t)_i^j - Q(N_r, 2N_t)_{i}^{2j}\|_{\infty} $$

If $q_i^j = q(r_i, t_j) = (p(r_i, t_j), s(r_i, t_j))^\top$ is the exact solution on the same mesh as $Q(N_r, N_t)$ then in the asymptotic limit, we anticipate that

$$ \max_{i,j} \|Q(N_r, N_t)_i^j - q_i^j\|_{\infty} = \mathcal{O}(\Delta r^2) + \mathcal{O}(\Delta t) $$

and we make the ansatz that

$$ Q(N_r, N_t)_i^j - q_i^j = H(r_i, t_j)/N_r^2 + K(r_i, t_j)/N_t $$

for $O(1)$ error functions $H$ and $K$. From this, we have

$$ Q(N_r, N_t)_i^j - Q(2N_r, N_t)_{2i}^j = \frac{3H(r_i, t_j)}{4N_r^2} $$

$$ Q(N_r, N_t)_i^j - Q(N_r, 2N_t)_{i}^{2j} = \frac{K(r_i, t_j)}{2N_t}. $$

Finally, assuming the meshes are refined enough, we have that

$$ \Lambda_r Q(N_r, N_t) = \frac{3H^*}{4N_r^2} $$

$$ \Lambda_t Q(N_r, N_t) = \frac{K^*}{2N_t} $$
Assuming there are $N \ge 2$ peaks, let them occur at the times $0 < t_1 < t_2 < \dots < t_N < t_\infty$. We define a proxy for the mean frequency as

$$ \hat{\omega} = \frac{N-1}{t_N - t_1}. $$

Then we define

$$ \omega_{\text{avg}} = \begin{cases} \hat{\omega}, & t_{\infty} \le t_N + \max_{1 \le i \le N-1} (t_{i+1} - t_i) \\ 0, & \text{otherwise.} \end{cases} \qquad (46) $$

Intuitively, if the system is continuing to oscillate, then the end of the computational time $t_\infty$, should occur before the estimated latest possible time for the $(N+1)$st peak (otherwise there should have been $N+1$ peaks). If this condition fails, we define the frequency to be 0, indicating the system stopped oscillating.

## A.4 Validation of Numerical Scheme

Throughout the validation and simulations, the initial conditions prescribed are those of Equations (30)-(32).

To validate the convergence rates, we shall use mesh refinements. Let $N_r$ and $N_t$ be the meshing parameters. Denote $Q(N_r, N_t)_i^j = (p_i^j, q_i^j)^\top$ for $0 \le i \le N_r - 1, 0 \le j \le N_t$ to be the numerical approximation at $r = r_i$ and $t = t_j$ over the chosen grid. We define the operations

$$ \Lambda_r Q(N_r, N_t) = \max_{\substack{0 \le i \le N_r - 1 \\ 0 \le j \le N_t}} \|Q(N_r, N_t)_i^j - Q(2N_r, N_t)_{2i}^j\|_{\infty} $$

$$ \Lambda_t Q(N_r, N_t) = \max_{\substack{0 \le i \le N_r - 1 \\ 0 \le j \le N_t}} \|Q(N_r, N_t)_i^j - Q(N_r, 2N_t)_{i}^{2j}\|_{\infty} $$

If $q_i^j = q(r_i, t_j) = (p(r_i, t_j), s(r_i, t_j))^\top$ is the exact solution on the same mesh as $Q(N_r, N_t)$ then in the asymptotic limit, we anticipate that

$$ \max_{i,j} \|Q(N_r, N_t)_i^j - q_i^j\|_{\infty} = \mathcal{O}(\Delta r^2) + \mathcal{O}(\Delta t) $$

and we make the ansatz that

$$ Q(N_r, N_t)_i^j - q_i^j = H(r_i, t_j)/N_r^2 + K(r_i, t_j)/N_t $$

for $O(1)$ error functions $H$ and $K$. From this, we have

$$ Q(N_r, N_t)_i^j - Q(2N_r, N_t)_{2i}^j = \frac{3H(r_i, t_j)}{4N_r^2} $$

$$ Q(N_r, N_t)_i^j - Q(N_r, 2N_t)_{i}^{2j} = \frac{K(r_i, t_j)}{2N_t}. $$

Finally, assuming the meshes are refined enough, we have that

$$ \Lambda_r Q(N_r, N_t) = \frac{3H^*}{4N_r^2} $$

$$ \Lambda_t Q(N_r, N_t) = \frac{K^*}{2N_t} $$
<table><caption>Table 8: Spatial convergence results.</caption><thead><tr><th>d</th><th>m<sub>r</sub></th><th>b<sub>r</sub></th><th>H*</th></tr></thead><tbody><tr><td>1</td><td>-1.987</td><td>6.275</td><td>708.292</td></tr><tr><td>2</td><td>-2.006</td><td>7.177</td><td>1744.776</td></tr><tr><td>3</td><td>-1.999</td><td>5.362</td><td>284.106</td></tr></tbody></table>

<table><caption>Table 9: Temporal convergence results</caption><thead><tr><th>d</th><th>m<sub>t</sub></th><th>b<sub>t</sub></th><th>K*</th></tr></thead><tbody><tr><td>1</td><td>-1.051</td><td>5.812</td><td>668.522</td></tr><tr><td>2</td><td>-0.999</td><td>4.282</td><td>144.800</td></tr><tr><td>3</td><td>-0.971</td><td>2.629</td><td>27.720</td></tr></tbody></table>

where $H^* = \sup_{[\rho, r_\infty] \times [0, t_\infty]} ||H||_\infty$ and $K^* = \sup_{[\rho, r_\infty] \times [0, t_\infty]} ||K||_\infty$.

In our validation, we focus our attention on simulations near the biologically relevant parameters of Table 5 and the phase diagrams of Figures 3-7. With $t_\infty = 80$, we choose $\gamma = 5$, $\sigma = 0.01$, $\delta = 0.05$, $\beta = 0.01$, $\eta = 1$, $\tau = 5$, and $\rho = 0.15$.

Numerically, we fix $N_t = 2^{14}$ and vary $N_r$ over the values $2^i$ where $i = 7, 8, 9, 10$ to compute $\Lambda_r$-values. Then we fix $N_r = 2^{10}$ and vary $N_t$ over the values $2^j$ where $j = 11, 12, 13, 14$ to compute $\Lambda_t$-values. We have

$$ \log \Lambda_r Q(N_r, N_t) = (\log(3/4) + \log(H^*)) - 2 \log N_r \qquad = b_r + m_r \log N_r $$

$$ \log \Lambda_t Q(N_r, N_t) = (\log(1/2) + \log(K^*)) - \log N_t \qquad = b_t + m_t \log N_t $$

for slopes $m_r$ and $m_t$ we expect to be close to $-2$ and $-1$, respectively, and intercepts $b_r$ and $b_t$. The values $H^*$ and $K^*$ can be found from $b_r$ and $b_t$. From lines of best fit, we estimate these values in Tables 8 and 9. Convergence plots are found in Figure 13.

These results suggest that to keep the asymptotic errors in both space and time below 0.005, we can use meshes as presented in Table 10.

# B Numerical Specifications

## B.1 Model Investigation

### B.1.1 Phase Diagram (Figures 3-7)

For a computational domain, we used $[\rho, r_\infty]$ with $\rho = 0.25$ and $r_\infty = 7.25$ for $r$, and $[0, t_\infty]$ with $t_\infty = 80$ for $t$. We used $r_* = (3\rho + r_\infty)/4$, $i_p = i_s = n_p = 1$. Our meshing parameters were those

<table><caption>Table 10: Estimates of mesh parameters to keep asymptotic errors $\le 0.005$.</caption><thead><tr><th>d</th><th>N<sub>t</sub></th><th>N<sub>r</sub></th></tr></thead><tbody><tr><td>1</td><td>150000</td><td>400</td></tr><tr><td>2</td><td>30000</td><td>600</td></tr><tr><td>3</td><td>6000</td><td>300</td></tr></tbody></table>
<table><caption>Table 8: Spatial convergence results.</caption><thead><tr><th>d</th><th>m<sub>r</sub></th><th>b<sub>r</sub></th><th>H*</th></tr></thead><tbody><tr><td>1</td><td>-1.987</td><td>6.275</td><td>708.292</td></tr><tr><td>2</td><td>-2.006</td><td>7.177</td><td>1744.776</td></tr><tr><td>3</td><td>-1.999</td><td>5.362</td><td>284.106</td></tr></tbody></table>

<table><caption>Table 9: Temporal convergence results</caption><thead><tr><th>d</th><th>m<sub>t</sub></th><th>b<sub>t</sub></th><th>K*</th></tr></thead><tbody><tr><td>1</td><td>-1.051</td><td>5.812</td><td>668.522</td></tr><tr><td>2</td><td>-0.999</td><td>4.282</td><td>144.800</td></tr><tr><td>3</td><td>-0.971</td><td>2.629</td><td>27.720</td></tr></tbody></table>

where $H^* = \sup_{[\rho, r_\infty] \times [0, t_\infty]} ||H||_\infty$ and $K^* = \sup_{[\rho, r_\infty] \times [0, t_\infty]} ||K||_\infty$.

In our validation, we focus our attention on simulations near the biologically relevant parameters of Table 5 and the phase diagrams of Figures 3-7. With $t_\infty = 80$, we choose $\gamma = 5$, $\sigma = 0.01$, $\delta = 0.05$, $\beta = 0.01$, $\eta = 1$, $\tau = 5$, and $\rho = 0.15$.

Numerically, we fix $N_t = 2^{14}$ and vary $N_r$ over the values $2^i$ where $i = 7, 8, 9, 10$ to compute $\Lambda_r$-values. Then we fix $N_r = 2^{10}$ and vary $N_t$ over the values $2^j$ where $j = 11, 12, 13, 14$ to compute $\Lambda_t$-values. We have

$$ \log \Lambda_r Q(N_r, N_t) = (\log(3/4) + \log(H^*)) - 2 \log N_r \qquad = b_r + m_r \log N_r $$

$$ \log \Lambda_t Q(N_r, N_t) = (\log(1/2) + \log(K^*)) - \log N_t \qquad = b_t + m_t \log N_t $$

for slopes $m_r$ and $m_t$ we expect to be close to $-2$ and $-1$, respectively, and intercepts $b_r$ and $b_t$. The values $H^*$ and $K^*$ can be found from $b_r$ and $b_t$. From lines of best fit, we estimate these values in Tables 8 and 9. Convergence plots are found in Figure 13.

These results suggest that to keep the asymptotic errors in both space and time below 0.005, we can use meshes as presented in Table 10.

# B Numerical Specifications

## B.1 Model Investigation

### B.1.1 Phase Diagram (Figures 3-7)

For a computational domain, we used $[\rho, r_\infty]$ with $\rho = 0.25$ and $r_\infty = 7.25$ for $r$, and $[0, t_\infty]$ with $t_\infty = 80$ for $t$. We used $r_* = (3\rho + r_\infty)/4$, $i_p = i_s = n_p = 1$. Our meshing parameters were those

<table><caption>Table 10: Estimates of mesh parameters to keep asymptotic errors $\le 0.005$.</caption><thead><tr><th>d</th><th>N<sub>t</sub></th><th>N<sub>r</sub></th></tr></thead><tbody><tr><td>1</td><td>150000</td><td>400</td></tr><tr><td>2</td><td>30000</td><td>600</td></tr><tr><td>3</td><td>6000</td><td>300</td></tr></tbody></table>
<figure><img src="image_48.png" /><figcaption>Figure 13: Plots of convergence with best fit.</figcaption></figure>
<figure><img src="image_46.png" /><figcaption>Figure 13: Plots of convergence with best fit.</figcaption></figure>
<table><caption>Table 11: Estimates of mesh parameters to keep asymptotic errors $\le 0.005$.</caption><thead><tr><th>d</th><th>N<sub>t</sub></th><th>N<sub>r</sub></th></tr></thead><tbody><tr><td>1</td><td>140000</td><td>400</td></tr><tr><td>2</td><td>120000</td><td>500</td></tr><tr><td>3</td><td>80000</td><td>300</td></tr></tbody></table>

presented in Table 10. The switch parameter $m = 10$ was used.

### B.1.2 Spatiotemporal Parameter Variations (Figures 8-10)

For a computational domain, we used $[\rho, r_\infty]$ (with $\rho$ as specified in the figures) with $r_\infty = \rho + 7$ for $r$, and $[0, t_\infty]$ with $t_\infty = 80$ for $t$. We used $r_* = (3\rho + r_\infty)/4$, $i_p = i_s = n_p = 1$. To ensure similar accuracy as in the phase diagram, we used meshing parameters outlined in Table 11. The switch parameter $m = 10$ was used.

### B.1.3 Comparison with Prior Work (Figure 11)

For a computational domain, we used $[0.25, r_\infty]$ with $r_\infty = \rho + 7$ for $r$, and $[0, t_\infty]$ with $t_\infty = 80$ for $t$. We used $r_* = (3\rho + r_\infty)/4$, $i_p = i_s = n_p = 1$. We used meshing parameters outlined in Table 10. The switch parameter $m = 10$ was used.

## B.2 Biological Investigation

### B.2.1 Generation of Tables 6-7 and Figure 12

For a computational domain, we used $[\rho, r_\infty]$ with $\rho$ found in Table 5 and $r_\infty = \rho + 7$ for $r$ and $[0, t_\infty]$ with $t_\infty = 80$ for $t$. We used $r_* = (3\rho + r_\infty)/4$. Our meshing parameters were those presented in Table 10. The switch parameter $m = 10$ was used. In this application, $i_p = i_s = 1$ again, but $n_p$ varied depending on the parameter that changed.

## References

[1] Goedert M, Spillantini MG. A century of Alzheimer's disease. science. 2006;314(5800):777-81.

[2] Davie CA. A review of Parkinson's disease. British medical bulletin. 2008;86(1):109-27.

[3] Akhtar A, Andleeb A, Waris TS, Bazzar M, Moradi AR, Awan NR, et al. Neurodegenerative diseases and effective drug delivery: A review of challenges and novel therapeutics. Journal of Controlled Release. 2021;330:1152-67.

[4] Kiaei M. New hopes and challenges for treatment of neurodegenerative disorders: Great opportunities for young neuroscientists. Basic and Clinical Neuroscience. 2013;4(1):3.

[5] Newby GA, Lindquist S. Blessings in disguise: biological benefits of prion-like mechanisms. Trends in cell biology. 2013;23(6):251-9.
<table><caption>Table 11: Estimates of mesh parameters to keep asymptotic errors ≤ 0.005.</caption><thead><tr><th>d</th><th>N<sub>t</sub></th><th>N<sub>r</sub></th></tr></thead><tbody><tr><td>1</td><td>140000</td><td>400</td></tr><tr><td>2</td><td>120000</td><td>500</td></tr><tr><td>3</td><td>80000</td><td>300</td></tr></tbody></table>

presented in Table 10. The switch parameter $m = 10$ was used.

### B.1.2 Spatiotemporal Parameter Variations (Figures 8-10)

For a computational domain, we used [$\rho, r_{\infty}$] (with $\rho$ as specified in the figures) with $r_{\infty} = \rho + 7$ for $r$, and [$0, t_{\infty}$] with $t_{\infty} = 80$ for $t$. We used $r_{*} = (3\rho + r_{\infty})/4$, $i_{p} = i_{s} = n_{p} = 1$. To ensure similar accuracy as in the phase diagram, we used meshing parameters outlined in Table 11. The switch parameter $m = 10$ was used.

### B.1.3 Comparison with Prior Work (Figure 11)

For a computational domain, we used [$0.25, r_{\infty}$] with $r_{\infty} = \rho + 7$ for $r$, and [$0, t_{\infty}$] with $t_{\infty} = 80$ for $t$. We used $r_{*} = (3\rho + r_{\infty})/4$, $i_{p} = i_{s} = n_{p} = 1$. We used meshing parameters outlined in Table 10. The switch parameter $m = 10$ was used.

## B.2 Biological Investigation

### B.2.1 Generation of Tables 6-7 and Figure 12

For a computational domain, we used [$\rho, r_{\infty}$] with $\rho$ found in Table 5 and $r_{\infty} = \rho + 7$ for $r$ and [$0, t_{\infty}$] with $t_{\infty} = 80$ for $t$. We used $r_{*} = (3\rho + r_{\infty})/4$. Our meshing parameters were those presented in Table 10. The switch parameter $m = 10$ was used. In this application, $i_{p} = i_{s} = 1$ again, but $n_{p}$ varied depending on the parameter that changed.

## References

[1] Goedert M, Spillantini MG. A century of Alzheimer's disease. science. 2006;314(5800):777-81.

[2] Davie CA. A review of Parkinson's disease. British medical bulletin. 2008;86(1):109-27.

[3] Akhtar A, Andleeb A, Waris TS, Bazzar M, Moradi AR, Awan NR, et al. Neurodegenerative diseases and effective drug delivery: A review of challenges and novel therapeutics. Journal of Controlled Release. 2021;330:1152-67.

[4] Kiaei M. New hopes and challenges for treatment of neurodegenerative disorders: Great opportunities for young neuroscientists. Basic and Clinical Neuroscience. 2013;4(1):3.

[5] Newby GA, Lindquist S. Blessings in disguise: biological benefits of prion-like mechanisms. Trends in cell biology. 2013;23(6):251-9.
[6] Ono K, Condron MM, Teplow DB. Structure–neurotoxicity relationships of amyloid β-protein oligomers. Proceedings of the National Academy of Sciences. 2009;106(35):14745-50.

[7] Tolar M, Abushakra S, Sabbagh M. The path forward in Alzheimer's disease therapeutics: Reevaluating the amyloid cascade hypothesis. Alzheimer's & Dementia. 2020;16(11):1553-60.

[8] Stefanis L. α-Synuclein in Parkinson's disease. Cold Spring Harbor perspectives in medicine. 2012;2(2):a009399.

[9] Prusiner S. Creutzfeldt-Jakob disease and scrapie prions. Alzheimer Disease & Associated Disorders. 1989;3(1):52-78.

[10] Hetz C, Zhang K, Kaufman RJ. Mechanisms, regulation and functions of the unfolded protein response. Nature reviews Molecular cell biology. 2020;21(8):421-38.

[11] Halliday M, Mallucci GR. Targeting the unfolded protein response in neurodegeneration: a new approach to therapy. Neuropharmacology. 2014;76:169-74.

[12] Fricker M, Tolkovsky AM, Borutaite V, Coleman M, Brown GC. Neuronal cell death. Physiological reviews. 2018;98(2):813-80.

[13] Apodaca J, Kim I, Rao H. Cellular tolerance of prion protein PrP in yeast involves proteolysis and the unfolded protein response. Biochemical and biophysical research communications. 2006;347(1):319-26.

[14] Hao W, Friedman A. Mathematical model on Alzheimer's disease. BMC systems biology. 2016;10(1):1-18.

[15] Puri IK, Li L. Mathematical modeling for the pathogenesis of Alzheimer's disease. PloS one. 2010;5(12):e15176.

[16] Lindstrom MR, Chavez MB, Gross-Sable EA, Hayden EY, Teplow DB. From reaction kinetics to dementia: A simple dimer model of Alzheimer's disease etiology. PLoS computational biology. 2021;17(7):e1009114.

[17] Desplats P, Lee HJ, Bae EJ, Patrick C, Rockenstein E, Crews L, et al. Inclusion formation and neuronal cell death through neuron-to-neuron transmission of α-synuclein. Proceedings of the National Academy of Sciences. 2009;106(31):13010-5.

[18] Bakshi S, Chelliah V, Chen C, van der Graaf PH. Mathematical biology models of Parkinson's disease. CPT: pharmacometrics & systems pharmacology. 2019;8(2):77-86.

[19] Pandya S, Zeighami Y, Freeze B, Dadar M, Collins DL, Dagher A, et al. Predictive model of spread of Parkinson's pathology using network diffusion. NeuroImage. 2019;192:178-94.

[20] Greer ML, Pujo-Menjouet L, Webb GF. A mathematical analysis of the dynamics of prion proliferation. Journal of theoretical biology. 2006;242(3):598-606.

[21] Salman S, Ahmed E. A mathematical model for Creutzfeldt Jacob Disease (CJD). Chaos, Solitons & Fractals. 2018;116:249-60.
[6] Ono K, Condron MM, Teplow DB. Structure–neurotoxicity relationships of amyloid β-protein oligomers. Proceedings of the National Academy of Sciences. 2009;106(35):14745-50.

[7] Tolar M, Abushakra S, Sabbagh M. The path forward in Alzheimer's disease therapeutics: Reevaluating the amyloid cascade hypothesis. Alzheimer's & Dementia. 2020;16(11):1553-60.

[8] Stefanis L. α-Synuclein in Parkinson's disease. Cold Spring Harbor perspectives in medicine. 2012;2(2):a009399.

[9] Prusiner S. Creutzfeldt-Jakob disease and scrapie prions. Alzheimer Disease & Associated Disorders. 1989;3(1):52-78.

[10] Hetz C, Zhang K, Kaufman RJ. Mechanisms, regulation and functions of the unfolded protein response. Nature reviews Molecular cell biology. 2020;21(8):421-38.

[11] Halliday M, Mallucci GR. Targeting the unfolded protein response in neurodegeneration: a new approach to therapy. Neuropharmacology. 2014;76:169-74.

[12] Fricker M, Tolkovsky AM, Borutaite V, Coleman M, Brown GC. Neuronal cell death. Physiological reviews. 2018;98(2):813-80.

[13] Apodaca J, Kim I, Rao H. Cellular tolerance of prion protein PrP in yeast involves proteolysis and the unfolded protein response. Biochemical and biophysical research communications. 2006;347(1):319-26.

[14] Hao W, Friedman A. Mathematical model on Alzheimer's disease. BMC systems biology. 2016;10(1):1-18.

[15] Puri IK, Li L. Mathematical modeling for the pathogenesis of Alzheimer's disease. PloS one. 2010;5(12):e15176.

[16] Lindstrom MR, Chavez MB, Gross-Sable EA, Hayden EY, Teplow DB. From reaction kinetics to dementia: A simple dimer model of Alzheimer's disease etiology. PLoS computational biology. 2021;17(7):e1009114.

[17] Desplats P, Lee HJ, Bae EJ, Patrick C, Rockenstein E, Crews L, et al. Inclusion formation and neuronal cell death through neuron-to-neuron transmission of α-synuclein. Proceedings of the National Academy of Sciences. 2009;106(31):13010-5.

[18] Bakshi S, Chelliah V, Chen C, van der Graaf PH. Mathematical biology models of Parkinson's disease. CPT: pharmacometrics & systems pharmacology. 2019;8(2):77-86.

[19] Pandya S, Zeighami Y, Freeze B, Dadar M, Collins DL, Dagher A, et al. Predictive model of spread of Parkinson's pathology using network diffusion. NeuroImage. 2019;192:178-94.

[20] Greer ML, Pujo-Menjouet L, Webb GF. A mathematical analysis of the dynamics of prion proliferation. Journal of theoretical biology. 2006;242(3):598-606.

[21] Salman S, Ahmed E. A mathematical model for Creutzfeldt Jacob Disease (CJD). Chaos, Solitons & Fractals. 2018;116:249-60.
[22] Adimy M, Babin L, Pujo-Menjouet L. Neuron Scale Modeling of Prion Production with the Unfolded Protein Response. SIAM Journal on Applied Dynamical Systems. 2022;21(4):2487-517.

[23] Adimy M, Chekroun A, Pujo-Menjouet L, Sensi M. A multigroup approach to delayed prion production. Discrete and Continuous Dynamical Systems - Series B. 2024;29.

[24] Trusina A, Tang C. The unfolded protein response and translation attenuation: a modelling approach. Diabetes, Obesity and Metabolism. 2010;12:27-31.

[25] Goriely A, Kuhl E, Bick C. Neuronal oscillations on evolving networks: dynamics, damage, degradation, decline, dementia, and death. Physical review letters. 2020;125(12):128102.

[26] Garzón DN, Castillo Y, Navas-Zuloaga MG, Darwin L, Hardin A, Culik N, et al. Dynamics of prion proliferation under combined treatment of pharmacological chaperones and interferons. Journal of Theoretical Biology. 2021;527:110797.

[27] Atkinson CJ, Zhang K, Munn AL, Wiegmans A, Wei MQ. Prion protein scrapie and the normal cellular prion protein. Prion. 2016;10(1):63-82.

[28] Selkoe DJ. Folding proteins in fatal ways. nature. 2003;426(6968):900-4.

[29] Shafiq M, Da Vela S, Amin L, Younas N, Harris DA, Zerr I, et al. The prion protein and its ligands: Insights into structure-function relationships. Biochimica et Biophysica Acta (BBA)-Molecular Cell Research. 2022;1869(6):119240.

[30] Rangel A, Race B, Striebel J, Chesebro B. Non-amyloid and amyloid prion protein deposits in prion-infected mice differ in blockage of interstitial brain fluid. Neuropathology and applied neurobiology. 2013;39(3):217-30.

[31] Moreno JA, Radford H, Peretti D, Steinert JR, Verity N, Martin MG, et al. Sustained translational repression by eIF2α-P mediates prion neurodegeneration. Nature. 2012;485(7399):507-11.

[32] Di Domenico F, Lanzillotta C. Chapter Two - The disturbance of protein synthesis/degradation homeostasis is a common trait of age-related neurodegenerative disorders. Advances in Protein Chemistry and Structural Biology. 2022;132:49-87.

[33] Martin JB. Applications of molecular genetics to restorative neurology. In: Principles and Practice of Restorative Neurology. Elsevier; 1992. p. 189-201.

[34] Bera K, Kiepas A, Godet I, Li Y, Mehta P, Ifemembi B, et al. Extracellular fluid viscosity enhances cell migration and cancer dissemination. Nature. 2022;611(7935):365-73.

[35] Young M, Carroad P, Bell R. Estimation of diffusion coefficients of proteins. Biotechnology and bioengineering. 1980;22(5):947-55.

[36] Hrabe J, Hrabětová S, Segeth K. A model of effective diffusion and tortuosity in the extracellular space of the brain. Biophysical journal. 2004;87(3):1606-17.

[37] Masel J, Jansen VA, Nowak MA. Quantifying the kinetic parameters of prion replication. Biophysical chemistry. 1999;77(2-3):139-52.
[22] Adimy M, Babin L, Pujo-Menjouet L. Neuron Scale Modeling of Prion Production with the Unfolded Protein Response. SIAM Journal on Applied Dynamical Systems. 2022;21(4):2487-517.

[23] Adimy M, Chekroun A, Pujo-Menjouet L, Sensi M. A multigroup approach to delayed prion production. Discrete and Continuous Dynamical Systems - Series B. 2024;29.

[24] Trusina A, Tang C. The unfolded protein response and translation attenuation: a modelling approach. Diabetes, Obesity and Metabolism. 2010;12:27-31.

[25] Goriely A, Kuhl E, Bick C. Neuronal oscillations on evolving networks: dynamics, damage, degradation, decline, dementia, and death. Physical review letters. 2020;125(12):128102.

[26] Garzón DN, Castillo Y, Navas-Zuloaga MG, Darwin L, Hardin A, Culik N, et al. Dynamics of prion proliferation under combined treatment of pharmacological chaperones and interferons. Journal of Theoretical Biology. 2021;527:110797.

[27] Atkinson CJ, Zhang K, Munn AL, Wiegmans A, Wei MQ. Prion protein scrapie and the normal cellular prion protein. Prion. 2016;10(1):63-82.

[28] Selkoe DJ. Folding proteins in fatal ways. nature. 2003;426(6968):900-4.

[29] Shafiq M, Da Vela S, Amin L, Younas N, Harris DA, Zerr I, et al. The prion protein and its ligands: Insights into structure-function relationships. Biochimica et Biophysica Acta (BBA)-Molecular Cell Research. 2022;1869(6):119240.

[30] Rangel A, Race B, Striebel J, Chesebro B. Non-amyloid and amyloid prion protein deposits in prion-infected mice differ in blockage of interstitial brain fluid. Neuropathology and applied neurobiology. 2013;39(3):217-30.

[31] Moreno JA, Radford H, Peretti D, Steinert JR, Verity N, Martin MG, et al. Sustained translational repression by eIF2α-P mediates prion neurodegeneration. Nature. 2012;485(7399):507-11.

[32] Di Domenico F, Lanzillotta C. Chapter Two - The disturbance of protein synthesis/degradation homeostasis is a common trait of age-related neurodegenerative disorders. Advances in Protein Chemistry and Structural Biology. 2022;132:49-87.

[33] Martin JB. Applications of molecular genetics to restorative neurology. In: Principles and Practice of Restorative Neurology. Elsevier; 1992. p. 189-201.

[34] Bera K, Kiepas A, Godet I, Li Y, Mehta P, Ifemembi B, et al. Extracellular fluid viscosity enhances cell migration and cancer dissemination. Nature. 2022;611(7935):365-73.

[35] Young M, Carroad P, Bell R. Estimation of diffusion coefficients of proteins. Biotechnology and bioengineering. 1980;22(5):947-55.

[36] Hrabe J, Hrabětová S, Segeth K. A model of effective diffusion and tortuosity in the extracellular space of the brain. Biophysical journal. 2004;87(3):1606-17.

[37] Masel J, Jansen VA, Nowak MA. Quantifying the kinetic parameters of prion replication. Biophysical chemistry. 1999;77(2-3):139-52.
[38] Beekes M, Baldauf E, Diringer H. Sequential appearance and accumulation of pathognomonic markers in the central nervous system of hamsters orally infected with scrapie. Journal of General Virology. 1996;77(8):1925-34.

[39] Barber TW, Brockway JA, Higgins LS. The density of tissues in and about the head. Acta neurologica scandinavica. 1970;46(1):85-92.

[40] Matthäus F. Diffusion versus network models as descriptions for the spread of prion diseases in the brain. Journal of theoretical biology. 2006;240(1):104-13.

[41] Peña J, Harris E. Dengue virus modulates the unfolded protein response in a time-dependent manner. Journal of Biological Chemistry. 2011;286(16):14226-36.

[42] Braitenberg V. Brain size and number of neurons: an exercise in synthetic neuroanatomy. Journal of computational neuroscience. 2001;10:71-7.

[43] Raine C. Characteristics of the Neuron; 1999. https://www.ncbi.nlm.nih.gov/books/NBK28209/.

[44] Stuart G, Spruston N, Häusser M. Dendrites. Oxford University Press; 2016.

[45] Rainov N, Tsuboi Y, Krolak-Salmon P, Vighetto A, Doh-Ura K. Experimental treatments for human transmissible spongiform encephalopathies: is there a role for pentosan polysulfate? Expert opinion on biological therapy. 2007;7(5):713-26.

[46] Wade J. Belfast man suffering from variant CJD passes away. The Journal. 2011. Available from: https://www.thejournal.ie/belfast-man-suffering-from-variant-cjd-passes-away-100653-Mar2011/.

[47] Yamasaki T, Suzuki A, Hasebe R, Horiuchi M. Comparison of the anti-prion mechanism of four different anti-prion compounds, anti-PrP monoclonal antibody 44B1, pentosan polysulfate, chlorpromazine, and U18666A, in prion-infected mouse neuroblastoma cells. PLoS One. 2014;9(9):e106516.

[48] Douglas Jr J. A survey of numerical methods for parabolic differential equations. In: Advances in computers. vol. 2. Elsevier; 1961. p. 1-54.
[38] Beekes M, Baldauf E, Diringer H. Sequential appearance and accumulation of pathognomonic markers in the central nervous system of hamsters orally infected with scrapie. Journal of General Virology. 1996;77(8):1925-34.

[39] Barber TW, Brockway JA, Higgins LS. The density of tissues in and about the head. Acta neurologica scandinavica. 1970;46(1):85-92.

[40] Matthäus F. Diffusion versus network models as descriptions for the spread of prion diseases in the brain. Journal of theoretical biology. 2006;240(1):104-13.

[41] Peña J, Harris E. Dengue virus modulates the unfolded protein response in a time-dependent manner. Journal of Biological Chemistry. 2011;286(16):14226-36.

[42] Braitenberg V. Brain size and number of neurons: an exercise in synthetic neuroanatomy. Journal of computational neuroscience. 2001;10:71-7.

[43] Raine C. Characteristics of the Neuron; 1999. https://www.ncbi.nlm.nih.gov/books/NBK28209/.

[44] Stuart G, Spruston N, Häusser M. Dendrites. Oxford University Press; 2016.

[45] Rainov N, Tsuboi Y, Krolak-Salmon P, Vighetto A, Doh-Ura K. Experimental treatments for human transmissible spongiform encephalopathies: is there a role for pentosan polysulfate? Expert opinion on biological therapy. 2007;7(5):713-26.

[46] Wade J. Belfast man suffering from variant CJD passes away. The Journal. 2011. Available from: https://www.thejournal.ie/belfast-man-suffering-from-variant-cjd-passes-away-100653-Mar2011/.

[47] Yamasaki T, Suzuki A, Hasebe R, Horiuchi M. Comparison of the anti-prion mechanism of four different anti-prion compounds, anti-PrP monoclonal antibody 44B1, pentosan polysulfate, chlorpromazine, and U18666A, in prion-infected mouse neuroblastoma cells. PLoS One. 2014;9(9):e106516.

[48] Douglas Jr J. A survey of numerical methods for parabolic differential equations. In: Advances in computers. vol. 2. Elsevier; 1961. p. 1-54.