September 25, 2017, Monday

# MBW:Extension to Mathematical Modeling of Alzheimer's Disease

The following is a replication of the work done by Ishwar K. Puri, Liwu Li in Mathematical Modeling for the Pathogenesis of Alzheimer's Disease. For more background information go to MBW: Mathematical Modeling of Alzheimer's Disease.

# Executive Summary

The following is a reproduction of the results of Puri et. al in Mathematical Modeling for the Pathogenesis of Alzheimer’s Disease [1]. Alzheimer's is one of the most prevalent causes of dementia and even death in people over 65 years old, however little is known about the development of the disease on a neurological level. The authors studied the mechanism of one of the proposed causes of Alzheimer's, the effect of the deposition of amyloid-$\beta$ particles forming senile plaques in the brain. Using a system of ordinary differential equations, it is possible to model one of the mechanisms taking place in the senile plaques, the activation of active microglia that in turn send out toxins and stimulate the death of the brain cells. The results from this model show that inflammatory activation of microglia is the dominant factor to the development of Alzheimer's and that amyloid-$\beta$ plays a principal role in the mechanism of activation of microglia, however is not itself the key player in the development of the disease. In terms of treatment, however, decreasing the presence of amyloid-$\beta$ is shown to decrease the cell death associated with Alzheimer's.

# Background

One of the most common causes of dementia is Alzheimer’s, a neurodegenerative disorder that is predominant in the senile population. It is estimated that 11% of people 65 years old and older and 32% of those 85 years old and older have Alzheimer’s. In addition to its predominance in society, this disorder has a high mortality rate, at the sixth leading cause of death in the United States. Just in 2012, more than \$216 billion were spent on unpaid healthcare for patients with Alzheimer's. [2] Clearly, there are a lot of societal costs. Even worse are the costs directly to patients and their families. Patients with Alzheimer's have tissue loss in the brain, resulting in brain shrinkage in a multitude of areas including the hipposcampus, which is responsible for storing long-term memories.

A key phenotype of patients that possess Alzheimer’s is the prevalence of amyloid-$\beta$ peptides deposited in the brain. It is believed that amyloid-$\beta$ peptides are present in humans of all ages, but as people age, the body loses ability to maintain an adequate level, and starts to deposit these proteins in senile plaques, or small masses of argyrophilic particles in the brain tissue. These plaques tend to be no more than 50 or 60 micrometers in diameter, but it appears that they still have a significant effect on brain activity. [3] These plaques are thought to stimulate the death of neurons, or brain cells that relay and store information.

However the mechanism of the development of this disorder has yet to be fully understood.Scientists have identified several key factors that are at play in the progression of Alzheimer's. One such factor that is known to have a strong contribution is genetic mutations. Currently, it is known that some of patients' Alzheimer's can be attributed to mutations in three specific genes. These genes are the ones responsible for the precursors of the amyloid-$\beta$ protein. Inheriting the mutations in those genes guarantees that a patient will develop Alzheimer's. However, less than 1% of all cases can be explained by this cause. The Apolipoprotein E-e4 (APOE-e4) Gene also appears to be prevalent in patients with Alzheimer's, but links have yet to be made as to the mechanism by witch it acts. [2]

It is proposed that the mechanism of development of Alzheimer's is centered around the presence of amyloid-$\beta$ particles, causing the inflammatory activation of specific neural cells, microglia. [1] These macrophages serve as immune cells in the brain, engulfing and breaking down foreign matter and debris, to ensure that the brain tissue is defended from infection. Microglia can be in one of two states, inflammatory or anti-inflammatory. When they are pro-inflammatory, or active, these cells secrete signaling molecules, reactive oxygen species, that may be toxic to neurons or other brain cells. [4] The pro-inflammatory state is coupled with a proliferating state of astrocytes that are responsible for maintaining neuron environments. When the astrocytes are in a quiescent state rather than the proliferating, the microglia display ani-inflammatory signals and do not send out any signals. More on the key players can be found at MBW: Mathematical Modeling of Alzheimer's Disease.

# Mathematical Model

Mechanism by Puri et al.
The proposed mechanism accounts for the feedback from amyloid-$\beta$, and the distinct states of microglia, astroglia, and neurons. The intercellular cross-talks are represented as the 16 pathways in the mechanism.

Two important and justified assumptions must be addressed:

1. Constant risk of neuronal death.

2. Spatiotemporal influence of diffusion is negligible.

However, within the scope of the small area of the senile plaques, both of these assumptions show external validity.

### Variable Definitions

1. ${\mathbf {N}}_{{s}}=$ Amount of neuron population that has survived the current time step.

2. ${\mathbf {N}}_{{d}}=$ Amount of neuron population that has died during the current time step.

3. ${\mathbf {A}}_{{q}}=$ Amount of astroglia population in a state of quiescence.

4. ${\mathbf {A}}_{{p}}=$ Amount of astroglia population in a state of proliferation.

5. ${\mathbf {M}}_{{2}}=$ Amount of normal microglia population in the resting anti-inflammatory state.

6. ${\mathbf {M}}_{{1}}=$ Amount of reactive microglia population in the active pro-inflammatory state.

7. ${\mathbf {A}}\beta =$ Number of amyloid-$\beta$ peptide molecules.

The initial conditions are give in Table 2 from Puri et al.

### Parameter Definitions

Table 1. Parameters
The proposed mechanism requires 16 rates for its description, and the model and analysis call for one additional rate. The required rates for the mechanism are given implicitly as $\alpha _{{i}}$ for $i=1,...,16$, which represent the cross-talks. The additional rate considered is $\alpha _{{r}}$, which represents the rate of amyloid-$\beta$ removal out of the system.

### Original System of Equations

The system is modeled with seven coupled ordinary differential equations. This system can also be represented as reaction rates or pathways between different states of cell populations:

This yields the following reaction rates:

(1.)

${\frac {d{\textbf {N}}_{s}}{dt}}=\alpha _{1}{\textbf {A}}_{q}-\alpha _{2}{\textbf {A}}_{p}-\alpha _{3}{\textbf {M}}_{1}$

(2.)

${\frac {d{\textbf {N}}_{d}}{dt}}=-{\frac {d{\textbf {N}}_{s}}{dt}}$

(3.)

${\frac {d{\textbf {A}}_{q}}{dt}}=\alpha _{4}{\textbf {M}}_{2}-\alpha _{5}{\textbf {M}}_{1}$

(4.)

${\frac {d{\textbf {A}}_{p}}{dt}}=-{\frac {d{\textbf {A}}_{q}}{dt}}$

(5.)

${\frac {d{\textbf {M}}_{2}}{dt}}=(\alpha _{6}+\alpha _{{11}}){\textbf {N}}_{s}-\alpha _{{10}}{\textbf {N}}_{d}+(\alpha _{7}+\alpha _{{12}}){\textbf {A}}_{q}-\alpha _{9}{\textbf {M}}_{1}+\alpha _{{14}}{\textbf {M}}_{2}-(\alpha _{8}+\alpha _{{13}}){\textbf {A}}\beta$

(6.)

${\frac {d{\textbf {M}}_{1}}{dt}}=-{\frac {d{\textbf {M}}_{2}}{dt}}$

(7a.)

${\frac {d{\textbf {A}}\beta }{dt}}=\alpha _{{15}}{\textbf {N}}_{s}-\alpha _{{16}}{\textbf {M}}_{2}$

These are identical to the original equations with the exception that one additional term is used to describe equation (7a.) as follows,

(7b.)

${\frac {d{\textbf {A}}\beta }{dt}}=\alpha _{{15}}{\textbf {N}}_{s}-\alpha _{{16}}{\textbf {M}}_{2}-\alpha _{r}{\textbf {A}}\beta$

Justification for the adapted model is presented in the Results section below.

# Sensitivity Analysis

Sensitivity analysis of the original system with respect to the parameter $\alpha _{1}$.

(1.1)

${\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {N}}_{s}}{dt}}\right)={\frac {\partial }{\partial \alpha _{1}}}\left(\alpha _{1}{\textbf {A}}_{q}-\alpha _{2}{\textbf {A}}_{p}-\alpha _{3}{\textbf {M}}_{1}\right)={\textbf {A}}_{q}+\alpha _{1}{\frac {\partial {\textbf {A}}_{q}}{\partial \alpha _{1}}}-\alpha _{2}{\frac {\partial {\textbf {A}}_{p}}{\partial \alpha _{1}}}-\alpha _{3}{\frac {\partial {\textbf {M}}_{1}}{\partial \alpha _{1}}}$

(2.1)

${\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {N}}_{d}}{dt}}\right)=-{\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {N}}_{s}}{dt}}\right)$

(3.1)

${\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {A}}_{q}}{dt}}\right)={\frac {\partial }{\partial \alpha _{1}}}\left(\alpha _{4}{\textbf {M}}_{2}-\alpha _{5}{\textbf {M}}_{1}\right)=\alpha _{4}{\frac {\partial {\textbf {M}}_{2}}{\partial \alpha _{1}}}-\alpha _{5}{\frac {\partial {\textbf {M}}_{1}}{\partial \alpha _{1}}}$

(4.1)

${\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {A}}_{p}}{dt}}\right)=-{\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {A}}_{q}}{dt}}\right)$

(5.1)

${\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {M}}_{2}}{dt}}\right)=$

${\frac {\partial }{\partial \alpha _{1}}}\left[(\alpha _{6}+\alpha _{{11}}){\textbf {N}}_{s}-\alpha _{{10}}{\textbf {N}}_{d}+(\alpha _{7}+\alpha _{{12}}){\textbf {A}}_{q}-\alpha _{9}{\textbf {M}}_{1}+\alpha _{{14}}{\textbf {M}}_{2}-(\alpha _{8}+\alpha _{{13}}){\textbf {A}}\beta \right]$

$=(\alpha _{6}+\alpha _{{11}}){\frac {\partial {\textbf {N}}_{s}}{\partial \alpha _{1}}}-\alpha _{{10}}{\frac {\partial {\textbf {N}}_{d}}{\partial \alpha _{1}}}+(\alpha _{7}+\alpha _{{12}}){\frac {\partial {\textbf {A}}_{q}}{\partial \alpha _{1}}}-\alpha _{9}{\frac {\partial {\textbf {M}}_{1}}{\partial \alpha _{1}}}+\alpha _{{14}}{\frac {\partial {\textbf {M}}_{2}}{\partial \alpha _{1}}}-(\alpha _{8}+\alpha _{{13}}){\frac {\partial {\textbf {A}}\beta }{\partial \alpha _{1}}}$

(6.1)

${\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {M}}_{1}}{dt}}\right)=-{\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {M}}_{2}}{dt}}\right)$

(7a.1)

${\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {A}}\beta }{dt}}\right)={\frac {\partial }{\partial \alpha _{1}}}\left(\alpha _{{15}}{\textbf {N}}_{s}-\alpha _{{16}}{\textbf {M}}_{2}\right)=\alpha _{{15}}{\frac {\partial {\textbf {N}}_{s}}{\partial \alpha _{1}}}-\alpha _{{16}}{\frac {\partial {\textbf {M}}_{2}}{\partial \alpha _{1}}}$

Then the corresponding adapted system is the same as above except that the last equation is defined as

(7b.1)

${\frac {\partial }{\partial \alpha _{1}}}\left({\frac {d{\textbf {A}}\beta }{dt}}\right)=\alpha _{{15}}{\frac {\partial {\textbf {N}}_{s}}{\partial \alpha _{1}}}-\alpha _{{16}}{\frac {\partial {\textbf {M}}_{2}}{\partial \alpha _{1}}}-\alpha _{r}{\frac {\partial {\textbf {A}}\beta }{\partial \alpha _{1}}}$.

# Results

## Cell Population Models

For the original model, an ordinary differential equation solver was used to produce the following plots in MatLab. The two plots in the first column are supposed to be identical to the Puri et al Figure 2 with $\alpha _{r}=1$, assuming that the results from Puri et al were indeed obtained by the model they described. The other two columns are also from the original model but with $\alpha _{r}=.1$ and $\alpha _{r}=.01$, respectively.

Notice that the amyloid-$\beta$ population is the same in each column despite varying the rate of amyloid-$\beta$ removal. This gave the motivation for the adapted model to include a term for amyloid-$\beta$ removal to equation (7).

Puri et al. Figure 2.

The following plots are from the adapted model. The two plots in the first column are identical to the Puri et al Figure 2, where the rate $\alpha _{r}=1$. The other two columns are also the adapted model but with $\alpha _{r}=.1$ and $\alpha _{r}=.01$, respectively.

Oddly enough, the Puri et al figure is unlike the MatLab version of the original system. However, the behavior of the original system makes sense by noticing that the first 6 equations represent one of two states for neurons, astroglia, and microglia but the seventh equation is only creating amyloid-$\beta$.

By the addition of a term for amyloid-$\beta$ removal in the seventh equation, the adapted model matches the Puri et al figure. We used a relaxed ODE solver in MatLab for stiff systems, yet the adapted model figure clearly has a negative cell population at about 19 years. We believe this is also true but less clear from the Puri et al figure.

## Sensitivity Comparison

The sensitivities for $\alpha _{1}$ of the original and adapted models calculated in MatLab are very different, using the equations described above. It is a little unclear how Puri et al preformed the sensitivity analysis. It seems they performed a finite difference, instead, by making small perturbations between $\pm 2.5\%$ in each $\alpha _{i}$ and determined after 20 simulation years. Even with the difference in calculation methods, the adapted model yields the same conclusions as cited in the paper by Puri et al. Below are the sensitivities calculated at simulation year 20.

Sensitivity of Neuron survival with respect to $\alpha _{1}$

As cited by Puri et al: $S(N_{s})=50000$

Original model MatLab calculations: $S(N_{s})=1.9811\times 10^{6}$

Adapted model MatLab calculations: $S(N_{s})=2.1206\times 10^{6}$

Sensitivity of Neuron death with respect to $\alpha _{1}$

As cited by Puri et al: $S(N_{d})=-50000$

Original model MatLab calculations: $S(N_{d})=-1.9811\times 10^{6}$

Adapted model MatLab calculations: $S(N_{d})=-2.1206\times 10^{6}$

Sensitivity of Microglia reactive with respect to $\alpha _{1}$

As cited by Puri et al: $S(M_{1})=-6000$

Original model MatLab calculations: $S(M_{1})=2.5213\times 10^{6}$

Adapted model MatLab calculations: $S(M_{1})=-2.8262\times 10^{5}$

Sensitivity of Microglia normal with respect to $\alpha _{1}$

As cited by Puri et al: $S(M_{2})=6000$

Original model MatLab calculations: $S(M_{2})=-2.5213\times 10^{6}$

Adapted model MatLab calculations: $S(M_{2})=2.8262\times 10^{5}$

Note that the exact value of the sensitivities is not as important as the sign and relative magnitude of the difference with respect to the other sensitivities. A positive sensitivity, $S(X_{j})$, implies that rate $\alpha _{i}$ corresponds to an increase of the $X_{j}$ population.

The following figure plots the sensitivities with respect to time for the original model. Notice these are incompatible with the behavior of the sensitivities cited by Puri et al, in both difference of magnitudes and sign. For example, the original sensitivities indicates that $\alpha _{1}$ promotes an increase of reactive microglia which is the opposite to that cited by Puri et al.

The following figure plots the sensitivities of the adapted model. Unlike the original model above, the behavior of the sensitivities match those cited by Puri et al.

It is possible to see that cell death and cell survival have opposite sensitivities, as expected, because the two states are in competition with each other, just as are the two states of the microglia.

## Variation in Parameters

Reactive Microglia with Respect to Changes in amyloid-$\beta$ Removal Rate
Proliferating Astroglia with Respect to Changes in amyloid-$\beta$ Removal Rate
Amyloid-beta with Respect to Changes in amyloid-$\beta$ Removal Rate

It is also possible to reproduce the findings of Puri et al. with respect to changes in each of the different populations when adjusting the values of different parameters. First, changes in the reactive microglia are studied when changing the rate of amyloid-$\beta$ removal. For higher values of amyloid-$\beta$ removal from the system, the amount of reactive microglia decreases. Looking at the proliferating astroglia, changes with respect to amyloid-$\beta$ removal as well, with increases in amyloid-$\beta$ removal resulting in the decrease of Astroglia proliferation. In addition, as expected, when the rate at which amyloid-$\beta$ is removed from the system increase, the amount of Amyloid-beta decreases.

This leads to the main results that, as currently suspected, removing amyloid-$\beta$ proteins from the brain would decrease cell death.

However, this only studies one of the mechanism at play. Looking at the rate of formation of amyloid-$\beta$ by the resting microglia, $\alpha _{{13}}$, it is clear that the amyloid-$\beta$ is just an intermediate actor. As the rate of creation of amyloid-$\beta$ increases, there is an increase in the production of the pro-inflammatory microglia:

Cell Death with Respect to Changes in amyloid-$\beta$ Formation Rate by Active Microglia
Amyloid-$\beta$ with Respect to Changes in amyloid-$\beta$ Formation Rate by Active Microglia

In turn, it is clear that as the rate of creation increases, the cell death increases and the reactive microglia population initially increases, but then decreases with time.

# Conclusions

The original model as described in the paper did not yield the same results as cited by Puri et al. However, the adapted model fits the results remarkably well. The figures from the adapted model are, to the best of our knowledge, identical to the figures by Puri et al. We believe our sensitivity calculations are more accurate than what was described by Puri et al., however both ways of calculating lead to the same conclusions, that the key factors at play in this mechanism are the active microglia and the amyloid-$\beta$ proteins. Yet again, the adapted model reproduces results closest to that of the paper.

Assuming the adapted model is correct and intended by the authors of the paper, our results support the conclusions. Namely, the importance of amyloid-$\beta$ and, more importantly, the role of reactive microglia as a major component in the pathogenesis of Alzheimer's disease. This model has shown that inflammatory activation of microglia is the dominant factor to the development of Alzheimer's and that amyloid-$\beta$ plays a principal role in the mechanism of activation of microglia, however is not itself the key player in the development of the disease. Even though the amyloid-$\beta$ removal is not physiologically valid, the authors used that to show that decreasing the presence of amyloid-$\beta$ is shown to decrease the cell death associated with Alzheimer's. Thus, a treatment targeting the amount of this protein or the rate of the formation of this protein by the anti-inflammatory microglia would be the best course of action for treatment of the disease.

# Code

### Main.m

%% Main:
% This program initiates the simulation and calls on Odesystem.m

clear all
clc
format long

%% Initial conditions:
% Declare global variables:
global y0 alpha;

% Preallocating vectors:
y0 = zeros(7,1);
alpha = zeros(17,1);

% Cell and molecule populations at time t(0) for some arbitrary volume:
y0(1) = 1e+4; % N_s    : Neuron survival.
y0(2) = 1e+2; % N_d    : Neuron death.
y0(3) = 1e+5; % A_q    : Astroglia quiescent.
y0(4) = 1e+3; % A_p    : Astroglia proliferation.
y0(5) = 1e+5; % M_2    : Microglia normal.
y0(6) = 1e+3; % M_1    : Microglia reactive.
y0(7) = 1e+3; % A_beta : Amyloid-beta molecules.

% Rates associated to the proposed mechanism. Units are [1/year]:
alpha(1) = 1e-5;
alpha(2) = 1e-3;
alpha(3) = 1e-2;
alpha(4) = 1e-4;
alpha(5) = 1e-2;
alpha(6) = 1e-2;
alpha(7) = 1e-4;
alpha(8) = 1e-2;
alpha(9) = 1e-2;
alpha(10) = 1e-2;
alpha(11) = 1e-2;
alpha(12) = 1e-4;
alpha(13) = 1e-2;
alpha(14) = 1e-4;
alpha(15) = 1;
alpha(16) = 1e-2;
alpha(17) = 1;

% I/O: If Toggle(1) is active then Toggle(2) must be commented out, and vice versa.
% Loop through three different Amyloid-beta removal rates for alpha_r:
%% Toggle(1): Begin

for j=1:3
if j==1
alpha(17) = 1; % alpha_r : Removal rate of Amyloid-beta from the system.
end

if j==2
alpha(17) = .1;
end

if j==3
alpha(17) = .01;
end

%% Toggle(1): End

% Loop through three different alpha_13 rates for impact of A_beta -> M_1:
%% Toggle(2): Begin

% for j=1:3
%     if j==1
%         alpha(13) = 1e-2; % alpha_13 :.
%     end
%
%     if j==2
%         alpha(13) = .1;
%     end
%
%     if j==3
%         alpha(13) = .5;
%     end

%% Toggle(2): End

%% Time:

% Move forward in time from 0 to 20. Units of time are years:
t0 = 0;
tf = 20;

%% Solve the system of difference equations with above initial conditions:

[t,y] = ode15s( @OdeSystem, [t0 tf], y0 );

%% Plots:
% The plots below are to be run with Toggle(1) active above unless
% stated otherwise (Plot 4 requires Toggle(2) active).

zero = t; % Straight line for comparison, y=0.

%% Plot 2:

% (2a.1) Plot for alpha_r = 1:
if j==1
%figure
subplot(2,3,1);
hold on
plot(t,y(:,1),'g'); % N_s : Neuron survival.
plot(t,y(:,6),'r'); % M_1 : Microglia reactive.
plot(t,y(:,7),'b'); % A_beta : Amyloid-beta molecules.
plot(t,zero,'-k', 'Linewidth',1);
title('\bf Populations with Amyloid-{\beta} removal: {\alpha_r = 1}','FontSize',16);
legend('N_s: Neuron Survival','M_1: Reactive Microglia','A{\beta}: Amyloid-\beta',2);
xlabel('\bf Time [years]','FontSize',12);
ylabel('\bf Populations of {N_s}, {M_1}, and {A{\beta}}','FontSize',12);
%axis([0 20 -100 11100]);
hold off
end

% (2a.2) Plot for alpha_r = .1:
if j==2
%figure
subplot(2,3,2);
hold on
plot(t,y(:,1),'g'); % N_s : Neuron survival.
plot(t,y(:,6),'r'); % M_1 : Microglia reactive.
plot(t,y(:,7),'b'); % A_beta : Amyloid-beta molecules.
plot(t,zero,'-k', 'Linewidth',1);
title('\bf {\alpha_r = .1}','FontSize',16);
%legend('N_s: Neuron Survival','M_1: Reactive Microglia','A{\beta}: Amyloid-\beta',2);
xlabel('\bf Time [years]','FontSize',12);
ylabel('\bf Populations of {N_s}, {M_1}, and {A{\beta}}','FontSize',12);
%axis([0 20 -100 15100]);
hold off
end

% (2a.3) Plot for alpha_r = .01:
if j==3
%figure
subplot(2,3,3);
hold on
plot(t,y(:,1),'g'); % N_s : Neuron survival.
plot(t,y(:,6),'r'); % M_1 : Microglia reactive.
plot(t,y(:,7),'b'); % A_beta : Amyloid-beta molecules.
plot(t,zero,'-k', 'Linewidth',1); % The line y = 0.
title('\bf {\alpha_r = .01}','FontSize',16);
%legend('N_s: Neuron Survival','M_1: Reactive Microglia','A{\beta}: Amyloid-\beta',2);
xlabel('\bf Time [years]','FontSize',12);
ylabel('\bf Populations of {N_s}, {M_1}, and {A{\beta}}','FontSize',12);
%axis([0 20 -100 15100]);
hold off
end

% (2b.1) Plot for alpha_r = 1:
if j==1
%figure
subplot(2,3,4);
hold on
plot(t,y(:,2),'g'); % N_d : Neuron death.
plot(t,y(:,4),'b'); % A_p : Astroglia proliferation.
title('\bf Neuron death and Astroglia proliferation:','FontSize',16);
legend('N_d: Neuron Death','A_p: Astroglia Proliferation',2);
xlabel('\bf Time [years]','FontSize',12);
ylabel('\bf Populations of {N_d}, {A_p}','FontSize',12);
hold off
end

% (2b.2) Plot for alpha_r = .1:
if j==2
%figure
subplot(2,3,5);
hold on
plot(t,y(:,2),'g'); % N_d : Neuron death.
plot(t,y(:,4),'b'); % A_p : Astroglia proliferation.
%title('\bf {\alpha_r = .1}','FontSize',16);
%legend('N_d: Neuron Death','A_p: Astroglia Proliferation',2);
xlabel('\bf Time [years]','FontSize',12);
ylabel('\bf Populations of {N_d}, {A_p}','FontSize',12);
hold off
end

% (2b.3) Plot for alpha_r = .01:
if j==3
%figure
subplot(2,3,6);
hold on
plot(t,y(:,2),'g'); % N_d : Neuron death.
plot(t,y(:,4),'b'); % A_p : Astroglia proliferation.
%title('\bf {\alpha_r = .01}','FontSize',16);
%legend('N_d: Neuron Death','A_p: Astroglia Proliferation',2);
xlabel('\bf Time [years]','FontSize',12);
ylabel('\bf Populations of {N_d}, {A_p}','FontSize',12);
hold off
end

%% Plot 3:
% Must comment out plots 1-2 above and 4-5 below to plot this section:
% Warning ~ Must plot 3a, 3b, 3c, and 3d ~>ONE AT A TIME<~ with All other
% plots commented out.

%% (3a) Reactive microglia (M_1) for three values of alpha_r:

%     hold on
%     if j==1
%
%         plot(t,y(:,6),'k');
%     end
%     if j==2
%         plot(t,y(:,6),'r');
%     end
%     if j==3
%         plot(t,y(:,6),'b');
%         title('Reactive Microglia for different rates of Amyloid-beta removal');
%         legend('{\alpha_r = 1}','{\alpha_r = .1}','{\alpha_r = .01}');
%         xlabel('Time [years]');
%         ylabel('Population of {M_1} for different {\alpha_r} values');
%         hold off
%     end

%% (3b) Astroglia proliferation (A_p) for three values of alpha_r:

%     hold on
%     for v=1:3
%
%         if j==1
%             plot(t,y(:,4),'k');
%         end
%         if j==2
%             plot(t,y(:,4),'r');
%         end
%         if j==3
%             plot(t,y(:,4),'b');
%         end
%
%     end
%     title('Astroglia proliferation for different rates of Amyloid-beta removal');
%     xlabel('Time [years]');
%     ylabel('Population of {A_p} for different {\alpha_r} values');
%     hold off

%% (3c) Amyloid-beta (A_beta) for three values of alpha_r:

%     hold on
%     for v=1:3
%
%         if j==1
%             plot(t,y(:,7),'k');
%         end
%         if j==2
%             plot(t,y(:,7),'r');
%         end
%         if j==3
%             plot(t,y(:,7),'b');
%         end
%
%     end
%     title('Amyloid-beta for different rates of Amyloid-beta removal');
%     xlabel('Time [years]');
%     ylabel('Population of {A\beta} for different {\alpha_r} values');
%     hold off

%% (3d) Neuron death (N_d) for three values of alpha_r:

%     hold on
%     for v=1:3
%
%         if j==1
%             plot(t,y(:,2),'k');
%         end
%         if j==2
%             plot(t,y(:,2),'r');
%         end
%         if j==3
%             plot(t,y(:,2),'b');
%         end
%
%     end
%     title('Neuron death for different rates of Amyloid-beta removal');
%     %legend('{\alpha_r = 1}','{\alpha_r = .1}','{\alpha_r = .01}');
%     xlabel('Time [years]');
%     ylabel('Population of {N_d} for different {\alpha_r} values');
%     hold off

%% Plot 4:
% Must comment out plots 1-3 above and 5 below to plot this section:
% Warning ~ Must plot 4a, 4b, and 4c one at a time with All other
% plots commented out.
% Must run these with Toggle(2) active above.

%% (4a) Plot of M_1 population for alpha_13 = .01:

%     hold on
%     if j==1
%         plot(t,y(:,6),'k','Linewidth',1.2);
%     end
%     if j==2
%         plot(t,y(:,6),'r','Linewidth',1.2);
%     end
%     if j==3
%         plot(t,y(:,6),'b','Linewidth',1.2);
%         title('Impact of Amyloid-beta on reactive Microglia for{M_1}');
%         legend('{M_1 } for {\alpha_{13} = .01}','{M_1 } for {\alpha_{13} = .1}','{M_1 } for {\alpha_{13} = .5}');
%         xlabel('Time [years]');
%         ylabel('Population of {M_1} for different {\alpha_{13}} values');
%         hold off
%     end

%% (4b) Plot of N_d population for alpha_13 = .1:

%     hold on
%     if j==1
%         plot(t,y(:,2),'k','Linewidth',1.2);
%     end
%     if j==2
%         plot(t,y(:,2),'r','Linewidth',1.2);
%     end
%     if j==3
%         plot(t,y(:,2),'b','Linewidth',1.2);
%         title('Impact of Amyloid-beta on reactive Microglia for{N_d}');
%         legend('{N_d } for {\alpha_{13} = .01}','{N_d } for {\alpha_{13} = .1}','{N_d } for {\alpha_{13} = .5}');
%         xlabel('Time [years]');
%         ylabel('Population of {N_d} for different {\alpha_{13}} values');
%         hold off
%     end

%% (4c) Plot of A_beta population for alpha_13 = .5:

%     hold on
%     if j==1
%         plot(t,y(:,7),'k','Linewidth',1.2);
%     end
%     if j==2
%         plot(t,y(:,7),'r','Linewidth',1.2);
%     end
%     if j==3
%         plot(t,y(:,7),'b','Linewidth',1.2);
%         title('Impact of Amyloid-beta on reactive Microglia for{A\beta}');
%         legend('{A\beta } for {\alpha_{13} = .01}','{A\beta } for {\alpha_{13} = .1}','{A\beta } for {\alpha_{13} = .5}');
%         xlabel('Time [years]');
%         ylabel('Population of {A\beta} for different {\alpha_{13}} values');
%         hold off
%     end

%% Plot 5:

% (1) Plot for alpha_r = 1:
%     if j==1
%         figure
%         hold on
%         plot(t,y(:,1),'--g');  % N_s : Neuron survival.
%         plot(t,y(:,2),'og');  % N_d : Neuron death.
%         plot(t,y(:,3),'--b');  % A_q : Astroglia quiescent.
%         plot(t,y(:,4),'+b');  % A_p : Astroglia proliferation.
%         plot(t,y(:,5),'--r'); % M_2 : Microglia normal.
%         plot(t,y(:,6),'+r');  % M_1 : Microglia reactive.
%         plot(t,y(:,7),'cy');  % A_beta : Amyloid-beta molecules.
%         title('Neuron survival and Neuron death with {\alpha_r = 1}');
%         %legend('N_s','N_d','A_q','A_p','M_2','M_1','A{\beta}',1);
%         xlabel('Time [years]');
%         ylabel('All Populations and Molecules');
%         hold off
%     end
%
%     % (2) Plot for alpha_r = .1:
%     if j==2
%         figure
%         hold on
%         plot(t,y(:,1),'--g');  % N_s : Neuron survival.
%         plot(t,y(:,2),'+g');  % N_d : Neuron death.
%         plot(t,y(:,3),'--b');  % A_q : Astroglia quiescent.
%         plot(t,y(:,4),'+b');  % A_p : Astroglia proliferation.
%         plot(t,y(:,5),'--r'); % M_2 : Microglia normal.
%         plot(t,y(:,6),'+r');  % M_1 : Microglia reactive.
%         plot(t,y(:,7),'cy');  % A_beta : Amyloid-beta molecules.
%         title('Neuron survival and Neuron death with {\alpha_r = .1}');
%         %legend('N_s','N_d','A_q','A_p','M_2','M_1','A{\beta}',1);
%         xlabel('Time [years]');
%         ylabel('All Populations and Molecules');
%         hold off
%     end
%
%     % (3) Plot for alpha_r = .1:
%     if j==3
%         figure
%         hold on
%         plot(t,y(:,1),'--g');  % N_s : Neuron survival.
%         plot(t,y(:,2),'+g');  % N_d : Neuron death.
%         plot(t,y(:,3),'--b');  % A_q : Astroglia quiescent.
%         plot(t,y(:,4),'+b');  % A_p : Astroglia proliferation.
%         plot(t,y(:,5),'--r'); % M_2 : Microglia normal.
%         plot(t,y(:,6),'+r');  % M_1 : Microglia reactive.
%         plot(t,y(:,7),'cy');  % A_beta : Amyloid-beta molecules.
%         title('Neuron survival and Neuron death with {\alpha_r = .01}');
%         %legend('N_s','N_d','A_q','A_p','M_2','M_1','A{\beta}',1);
%         xlabel('Time [years]');
%         ylabel('All Populations and Molecules');
%         hold off
%     end

end


### OdeSystem.m

function dydt = OdeSystem(t,y0)
global alpha

N_s = y0(1);    % Neuron survival         : N_s(0) = 1.0e+4;
N_d = y0(2);    % Neuron death            : N_d(0) = 1.0e+2;
A_q = y0(3);    % Astroglia quiescent     : A_q(0) = 1.0e+5;
A_p = y0(4);    % Astroglia proliferation : A_p(0) = 1.0e+3;
M_2 = y0(5);    % Microglia normal        : M_2(0) = 1.0e+5
M_1 = y0(6);    % Microglia reactive      : M_1(0) = 1.0e+3;
A_beta = y0(7); % Amyloid-beta molecules  : A_beta(0) = 1.0e+3;

% Preallocate the solution vector:
dydt = zeros(7,1);

% Definitions of the seven coupled rate equations:
dydt(1) = alpha(1)*A_q - alpha(2)*A_p - alpha(3)*M_1;
dydt(2) = -(dydt(1));

dydt(3) = alpha(4)*M_2 - alpha(5)*M_1;
dydt(4) = -(dydt(3));

dydt(5) = (alpha(6) + alpha(11))*N_s - alpha(10)*N_d + (alpha(7) + alpha(12))*A_q ...
- alpha(9)*M_1 + alpha(14)*M_2 - (alpha(8) + alpha(13))*A_beta;

dydt(6) = -(dydt(5));

%% I/O: Can only run Main with either (7a) active and (7b) commented out, or vice versa:

% (7a) Original model:
%dydt(7) = alpha(15)*N_s - alpha(16)*M_2;

dydt(7) = alpha(15)*N_s - alpha(16)*M_2 - alpha(17)*A_beta;

% Note: The adapted model has one additional term required in dydt(7) for the rate at
% which Amyloid-beta is removed from the system; ie: -alpha(17)*A_beta.

return


### SensitivityMain.m

%% Sensitivity for alpha_1:
% This program initiates the simulation and calls on SensitivitySystem.m

clear all
clc
format long

%% Initial conditions:
% Declare global variables:
global y0 alpha;

% Preallocating vectors:
y0 = zeros(7,1);
alpha = zeros(17,1);

% Sensitivities of cell populations at time t(0) w.r.t alpha_1:
y0(1) = 1e+5;  % Value of initial A_q population from derivation of new system.
y0(2) = -1e+5; % This is derived from above.
y0(3) = 0;
y0(4) = 0;
y0(5) = 0;
y0(6) = 0;
y0(7) = 0;

% Rates associated to the proposed mechanism. Units are [1/year]:
alpha(1) = 1e-5;
alpha(2) = 1e-3;
alpha(3) = 1e-2;
alpha(4) = 1e-4;
alpha(5) = 1e-2;
alpha(6) = 1e-2;
alpha(7) = 1e-4;
alpha(8) = 1e-2;
alpha(9) = 1e-2;
alpha(10) = 1e-2;
alpha(11) = 1e-2;
alpha(12) = 1e-4;
alpha(13) = 1e-2;
alpha(14) = 1e-4;
alpha(15) = 1;
alpha(16) = 1e-2;
alpha(17) = 1;

%% Time:

% Move forward in time from 0 to 20. Units of time are years:
t0 = 0;
tf = 20;

%% Solve the system of difference equations with above initial conditions:

[t,y] = ode15s( @SensitivitySystem, [t0 tf], y0 );

figure
hold on
plot(t,y(:,1),'--g');  % dNs : d(Neuron survival).
plot(t,y(:,2),'og');  % dNd : d(Neuron death).
plot(t,y(:,5),'--r'); % dM2 : d(Microglia normal).
plot(t,y(:,6),'+r');  % dM1 : d(Microglia reactive).
title('\bf Sensitivity of cell populations w.r.t.  {\alpha_1}','FontSize',16);
legend('S(Ns)','S(Nd)','S(M2)','S(M1)',1);
xlabel('\bf Time [years]','FontSize',12);
ylabel('\bf Sensitivity to {\alpha_1}','FontSize',12);
hold off

%% Display S(Ns), S(Nd), S(M2), S(M1) with respect to alpha_1 at year 20:

[n m] = size(y);

disp('Sensitivity of Neuron survival with respect to alpha_1:');
SNs = y(n,1)
disp('     ');

disp('Sensitivity of Neuron death with respect to alpha_1:');
SNd = y(n,2)
disp('     ');

disp('Sensitivity of Microglia reactive with respect to alpha_1:');
SM1 = y(n,6)
disp('     ');

disp('Sensitivity of Microglia normal with respect to alpha_1:');
SM2 = y(n,5)
disp('     ');


### SensitivitySystem.m

function dydt = SensitivitySystem(t,y0)
global alpha

% This is for the population sensitivities with respect to
% alpha_1
dNs = y0(1);    % d(Neuron survival)         : dNs = 1.0e+5;
dNd = y0(2);    % d(Neuron death)            : dNd = -1.0e+5;
dAq = y0(3);    % d(Astroglia quiescent)     : dAq = 0;
dAp = y0(4);    % d(Astroglia proliferation) : dAp = 0;
dM2 = y0(5);    % d(Microglia normal)        : dM2 = 0;
dM1 = y0(6);    % d(Microglia reactive)      : dM1 = 0;
dAb = y0(7);    % d(Amyloid-beta molecules)  : dAb = 0;

Aq = 1e+5;

% Preallocate the solution vector:
dydt = zeros(7,1);

% Definitions of the seven coupled rate equations using
% d/d(alpha_1)(d(Population)/dt) :
dydt(1) = Aq + alpha(1)*dAq - alpha(2)*dAp - alpha(3)*dM1;
dydt(2) = -(dydt(1));
dydt(3) = alpha(4)*dM2 - alpha(5)*dM1;
dydt(4) = -(dydt(3));

dydt(5) = (alpha(6) + alpha(11))*dNs - alpha(10)*dNd + (alpha(7) + alpha(12))*dAq ...
- alpha(9)*dM1 + alpha(14)*dM2 - (alpha(8) + alpha(13))*dAb;

dydt(6) = -(dydt(5));

%% I/O: Can only run SensitivityMain with either (7a) active and (7b) commented out, or vice versa:

% (7a) Original model:
%dydt(7) = alpha(15)*dNs - alpha(16)*dM2;

dydt(7) = alpha(15)*dNs - alpha(16)*dM2 - alpha(17)*dAb;

% Note: one additional term was required in dydt(7) for the rate at
% which Amyloid-beta is removed from the system; ie: -alpha(17)*A_beta.

return


# References

1. Puri IK, Li L. (2010) Mathematical Modeling for the Pathogenesis of Alzheimer’s Disease. PLoS ONE 5(12): e15176. doi:10.1371/journal.pone.0015176
2. "2013 Alzheimer’s Disease Facts and Figures." Alzheimer’s & Dementia 9.2 (2013): 246-249. http://dx.doi.org/10.1016/j.jalz.2013.02.003.
3. Roth, Martin, Bernard E. Tomlinson, and Gary Blessed. "Correlation between scores for dementia and counts of ‘senile plaques’ in cerebral grey matter of elderly subjects." (1966): 109-110.
4. Smith, Joshua A., Arabinda Das, Swapan K. Ray, and Naren L. Banik. "Role of pro-inflammatory cytokines released from microglia in neurodegenerative diseases." Brain research bulletin 87.1 (2012): 10-20.