Mathematical model on Alzheimer's disease.

BACKGROUND: Alzheimer disease (AD) is a progressive neurodegenerative disease that destroys memory and cognitive skills. AD is characterized by the presence of two types of neuropathological hallmarks: extracellular plaques consisting of amyloid β-peptides and intracellular neurofibrillary tangles of hyperphosphorylated tau proteins. The disease affects 5 million people in the United States and 44 million world-wide. Currently there is no drug that can cure, stop or even slow the progression of the disease. If no cure is found, by 2050 the number of alzheimer's patients in the U.S. will reach 15 million and the cost of caring for them will exceed $ 1 trillion annually.
RESULTS: The present paper develops a mathematical model of AD that includes neurons, astrocytes, microglias and peripheral macrophages, as well as amyloid β aggregation and hyperphosphorylated tau proteins. The model is represented by a system of partial differential equations. The model is used to simulate the effect of drugs that either failed in clinical trials, or are currently in clinical trials.
CONCLUSIONS: Based on these simulations it is suggested that combined therapy with TNF- α inhibitor and anti amyloid β could yield significant efficacy in slowing the progression of AD.

Background

AD is the most common form of dementia. The disease is an irreversible, progressive, brain disorder that destroys memory and cognitive skills, and eventually the ability to carry out even the simplest tasks. While the genetic inheritability of AD is in the range of 50 –80% [1, 2], the cause of the disease is mostly unknown. The disease strikes ageing people typically 65 or older, and twice more women than men. In 2015 there were more than 5 million people in the United States with AD, and 44 millions world-wide [3]. The cost of caring for AD patients in the U.S. was estimated at $226 billions for 2015 [3].

AD is characterized by the presence of two types of neuropathological hallmarks: extracellular plaques and intracellular neurofibrillary tangles (NFTs). The extracellular plaques consist primarily of amyloid β-peptide (A β) deposits. The NFTs are intraneural aggregation of hyperphosphorylated tau proteins. Reactive oxygen species (ROS) appears to be one of the early events in the progression of the disease [4]. Amyloid precursor protein (APP) on neurons membrane constitutively shed A β peptides [5]. High levels of ROS promote abnormal deposition of A β [4, 6]. Tau protein in the central nervous system (CNS) is predominantly expressed in neurons; its main role is to promote microtubles assembly and stability. Glycogen synthase kinase-type 3(GSK-3) is activated by the abnormally produced A β, and it mediates the hyperphosphorylation of tau proteins [4, 6–9].

The hyperphosphorylated tau proteins cause microtuble depolymerization and destruction, as they aggregate to form neurofibrillary tangles. This results in neuronal death and release of the NFTs to the extracellular environment [4, 10].

The non-neuronal cells in the brain consist of cells that support neurons directly, mostly astrocytes, and immune cells.

Microglias are the resident macrophages in the brain. They constitute the main active immune cells in the brain. They are activated by soluble A β oligomers which build up from the A β deposits [11, 12].

Astrocytes are in close proximity to neurons. They support neuronal cross-talk, and mediate the transport of nutrients from the blood to neurons. Astrocytes are activated primarily by TNF- α, but also by A β [10, 13–16]. Activated astrocytes produce A β, but at a smaller rate than neurons [16]. Activated astrocytes also produce MCP-1, which attracts monocytes from the blood into the plaques [17-19]. The monocytes differentiate into proinflammatory macrophages, , but may then change phenotype into anti-inflammatory macrophage. Activated microglias have two phenotypes: proinflammatory M 1 macroglia and anti-inflammatory M 2 macroglia [12, 20]. Macrophages have a major role in A β clearance [12, 20], but activated microglia are poorly phagocytic for A β compared to peripheral macrophages [21]. M 1 and macrophages are neurotoxic; they produce proinflammatory cytokines TNF- α, IL-6, IL-12 and IL-1 β [20, 22, 23]. M 2 microglias and peripheral macrophages produce anti-inflammatory cytokines IL-10, IL-13, IL-4 and TGF- β [20]. The neuronal stress caused by the proinflammatory cytokines, is resisted by IL-10, IL-13 and IL-4, but nevertheless it contributes to neuronal damage and death [20, 22, 23].

There are currently no drugs that can cure AD, or stop its progression. Many clinical trials of drugs aimed at preventing or clearing the A β and tau pathology have failed to demonstrate efficacy [24-27]. Currently the only treatment of AD is by medications that are used to treat the symptoms of the disease.

The role of TGF- β is somewhat controversial [28]. On one hand, TGF- β provides protection against neuroninflammation and neurondegeneration [29-34], but on the other hand, TGF- β-induced TIAF1 interacts with amyloid fibrils to favorably support plaque formation [28], and blocking TGF- β-smad2/3 in peripheral macrophages mitigates AD pathology [35].

Figure 1 is a schematics of the network associated with the progression of AD. Figure 1 a shows the network within a neuron which leads from ROS to NFTs and the destruction of microtubules. Figure 1 b shows the network of activated cells, microglia, astrocyte and monocyte-derived macrophages and their effect on neurons and their microenvironment.

Fig. 1

Schematic network in AD: a Amyloid precurser protein (APP) sheds Amyloid β peptides. ROS promotes abnormal production of A β [5, 6], which activates GSK-3 [4, 6, 8]. Activated GSK-3 mediates hyperphosphrylation of tau proteins [4, 6], which results in the formation of NFTs [10] and destruction of microtubules [4, 10], leading to neuron death. b Astrocytes are activated by [10, 16] and TNF- α [14, 15], and they produce MCP-1 [17–19], which attracts macrophages into the tissue [17, 19]. NFT activates microglias [10, 13, 15]. Activated proinflammatory microglias and microphages produce TNF- α and other proinflammatory cytokines [20, 22, 23], while anti-inflammatory microglias and macrophages produce IL-10 and other anti-inflammatory cytokines [20, 22, 23]. Dead neurons release A β and NFTs, and soluble A β oligomers activate microglia [11, 12]. Activated astrocytes secrete A β [16]. A β deposit is reduced through endocytosis by microglia and macrophages [12, 20]

In this paper we develop a mathematical model of AD. The model is represented by a system of partial differential equations (PDEs) based on Fig. 1. For simplicity we represent all the proinflammatory cytokines by TNF- α, and all the anti-inflammatory cytokines by IL-10.

We shall use our model to conduct in silico trials with several drugs: TNF- α inhibitor, anti-A β drug, MCP-1 inhibitor, and injection of TGF- β. Simulations of the model show that continuous treatment with TNF- α inhibitor yields a slight decrease the death of neurons, and anti-A β drug yields a slight decrease in the aggregation of A β over 10 years period, while the benefits from injection of TGF- β and MCP-1 inhibitor drugs are negligible. This suggests that clinical trials consider combination therapy with TNF- α and anti-A β drugs.

We note that Fig. 1 does not display neurites: the projections of axons and dendrites from the body of neurons. It is known that the aggregations of A β mediate rapid disruption of synaptic plasticity and memory [36-39]. Thus the progression of AD in terms of reduction in dendritic complexity and synaptic dysfunction will not be considered in the present paper.

We conclude the Introduction by mentioning earlier mathematical models which deal with some aspects of AD: A β polymerization [40], A β plaque formation and the role of prions interacting with A β [41, 42], linear cross-talk among brain cells and A β [43], and the influence of SORLA on AD progression [44, 45].

Methods

Mathematical model

Model’s variables

The mathematical model is based on Fig. 1 and is represented by a system of partial differential equations. Table 1 lists the variables used in the model.

Table 1

The variables of the model; concentration and densities are in units of g/c m 3 for cells and g/m l for cytokines

ROS (R):	Reactive oxygen species	GSK-3 (G):	Glycogen synthase kinase-type 3
Aβi:	Amyloid β inside neurons	Aβo:	Amyloid β outside neurons
NFT (F i):	Neuronfibrillary tangle inside neurons	NFT (F o):	Neuronfibrillary tangle outside neurons
APP (A P):	Amyloid precursor protein	A βO (A O):	Amyloid β oligomer (soluble)
TNF- α (T α):	Tumor necrosis factor alpha	TGF- β (T β):	Transforming growth factor beta
IL-10 (I 10):	Interleukin 10	P:	MCP-1
M 1:	Proinflammatory microglias	M 2:	Anti-inflammatory microglias
MG (M G):	Microglias	N:	Live neurons
A :	Astrocytes	N d:	Dead neurons
M^_1	Peripheral proinflammatory macrophages	M^_2:	Peripheral anti-inflammatory macrophages
τ	hyperphosphorylated tau protein	H	High mobility group box 1 (HMGB1)

Equations for A β

The amyloid- β within neurons, , are constitutively released from APP at a rate and are degraded at a rate . Under reactive oxidative stress, R, is overproduced. Hence the equation for is given by

where N 0 is the reference density of the neuron cells in the brain.

The extracellular amyloid- β peptides satisfy the following equation:

where is a Michaelis-Menten coefficient. Neurons die at a rate , thereby releasing their . Hence they contribute to the growth rate of , which is the first term on the right-hand side of Eq. (2). The second term on the right-hand side of Eq. (2) represents A β constitutively released from APP [5], and the third term accounts for A β released by activated astrocytes [16]; A 0 is the reference density of the astrocyte cells in the brain. is cleared primarily by peripheral macrophages and , but also by activated microglias M 1 and M 2, so [21], and M 1 are more effective in clearing than and M 2 [46, 47] so 0≤θ<1. APP on live neurons shed A β peptides both inside the neurons (as ) and outside the neurons (as ). We assume that most are produced from dead neurons. Hence, in Eq. (2), we neglected the production of by live neurons. We also assumed that ROS increases primarily the A β that are within live neurons, and thus neglected the increase of by ROS.

Equation for

Tau protein is constitutively produced at some rate λ . We assume that when production exceeds a threshold , GSK-3 becomes activated and it mediates hyperphosphorylation of tau. In steady state, the difference is proportional to R. Hence the equation for tau is given by:

We assume that initially we already have a disease state. Thus, in particular, the tau proteins are already hyperphophorylated and ROS induces increases in the production of these proteins.

Equations for NFT

The NFTs in neurons (F ) are formed from the hyperphosphorylated tau proteins [4, 6–9], and they are released to the extraceullar space (and are then labeled F 0) when the neurons die [4, 10]. Hence,

Equation for neurons

Hyberphosphorated tau proteins, forming neurofibrillary tangles, cause microtubles depolymerization and destruction, resulting in neuron death [4, 6–9]. Neuron death is also caused by stress from proinflammatory cytokines which is, however, resisted by anti-inflammatory cytokines [20, 22, 23]. For simplicity we represent all the proinflammatory cytokines by TNF- α and all the anti-inflammatory cytokines by IL-10. Hence the equation for N takes the following form:

where the death rates of N caused by F and T are assumed to depend on their saturation levels.

Equation for astrocytes

Astrocytes are activated primarily by extracellular TNF- α [14, 15], but also by [10, 16], so that

Equation for dead neurons

The equation for dead neurons, N , is given by

where is a Michaelis-Menten coefficient. The first two terms on the right-hand side arise from the death of N cells. The last two terms account for the clearance of N by microglias and peripheral macrophages [48].

Equation for A O

The A βO are soluble A β oligomers and they can diffuse throughout the brain tissue [49, 50]. Their density A satisfies the equation:

where is the rate by which the A are formed from the extracellular amyloid β peptides, and accounts for the diffusion of A .

Equation for HMGB-1

In general, when cell death occurs through necrosis, dying cells release HMGB-1 [51]. In AD, HMGB-1 is produced by dying neurons [52-54]. Hence,

Equations for activated microglias

Activated microglias have two phenotypes: proinflammatory M 1 and anti-inflammatory M 2. They satisfy the following equations:

where and .

Microglias can travel in the brain [55]. Activated microglias are chemoattracted to dead neurons [10, 13, 15], more precisely, to the cytokines HMGB-1 produced by N , and this is represented by the second term of the left-hand side of Eqs. (11), (12). Microglias are activated by extracelluar NFTs [10, 13, 15], and by soluble oligomers A [11, 12]. They become of M 1 phenotype under proinflammatory signals from TNF- α, and of M 2 phenotype under anti-inflammatory signals from IL-10. These facts are expressed by the first term on the right-hand sides of Eqs. (11), (12); is the ratio by which the activated microglias become M 1 macrophages, and is the ratio by which activated microglias become M 2 macrophages. The parameter β reflects the ratio of proinflammatory/anti-inflammatory environment, as determined by the relative ‘strength’ of T v.s. I 10.

In addition, there is a transition M 1→M 2 under the TGF- β signaling [32], which is accounted by the second term on the right-hand side of these equations.

Equations for macrophages

Peripheral macrophages are differentiated from monocytes which migrate through the blood vessels. They satisfy a flux condition on the boundary of the blood vessels, where n is the outward normal, M 0 is the density of the monocytes in the brain capillaries, and is a function which depends on the concentration of MCP-1. By averaging these fluxes from blood vessels, we can represent (as in [56]) the immigration of macrophages into the brain tissue by a term . We assume that the incoming macrophages divide into and phenotype depending on the relative concentrations of TNF- α and IL-10 [47]. Macrophages can also change phenotype to macrophages under signaling by TGF- β. We finally note that because of the blood-brain barrier (BBB) we do not include diffusion of peripheral macrophages, but we do include chemotaxis by amyloid- β plaques or, more specifically, by the soluble A [17-19]. Hence peripheral macrophages satisfy the following equations:

where and [56].

Equations for TGF- , TNF- MCP-1 and IL-10

T and IL-10 are produced by M 2 microglia and macrophages. TNF- α is produced by proinflammatory macrophages M 1 and . Hence the equations for T , T and I 10 have the following form:

MCP-1 is produced by activated astrocytes [17-19] and by microglias [12], which are assumed to be of M 2 phenotype. Hence

The estimates of parameters in Eqs. (1)–(18) are given in “Appendix".

Results and discussions

We simulate the model (1)-(18) in a rectangular domain Ω={(x,y),0≤x≤1,0≤y≤1}. We assume that

We take initial values, for each variable X, to be below (or above) the expected steady state for X, if X is expected to grow (or decrease) with the progression of the disease. A specific choice is given below, but the simulations do not change, after a short time, with other choices:

We also prescribe the value of ROS in Eqs. (1), (3) by

Figure 2 shows the average density of all the 18 variables of the model over a period of 10 years. We first observe that, for all the species that tend to a steady state in Fig. 2, the steady states are approximately the same as those that we assumed in estimating some of the model parameters. Thus the steady state values of and I 10 are approximately equal to the values assumed in “Appendix”. We conclude that estimates of the parameters which were based on steady state assumptions on macrophages, microglias and the half-saturation parameters are consistent with the simulation results.

Fig. 2

Average concentration of cytokines and average density of cells. All the parameters are as in Tables 2 and 3

Table 2

Parameters’ description and value

Parameter	Description	Value
D_A_O	Diffusion coefficient of A βO	4.32×10−2 c m 2 day −1 estimated
D H	Diffusion coefficient of HMGB-1	8.11×10−2 c m 2 day −1 estimated
D_T_α	Diffusion coefficient for TNF- α	6.55×10−2 c m 2 day −1 estimated
d_T_β	Diffusion coefficient of TGF- β	6.55×10−2 c m 2 day −1 estimated
D_I_10	Diffusion coefficient of IL-10	6.04×10−2 c m 2 day −1 estimated
D P	Diffusion coefficient of MCP-1	1.2×10−1 c m 2 day −1 estimated
λβi	Production rate of Aβi	9.51×10−6 g/ml/day estimated
λ N	Production rate of Aβo by neuron	8×10−9 g/ml/day estimated
λ A	Production rate of Aβo by astrocytes	8×10−10 g/ml/day estimated
λ τ0	Production rate of tau proteins in health	8.1×10−11 g/ml/day estimated
λ τ	Production rate of tau proteins by ROS	1.35×10−11 g/ml estimated
λ F	Production rate of NFT by tau	1.662×10−3/day estimated
λ_AT_α	Production/activation rate of astrocytes by TNF- α	1.54/day estimated
λ_AAβo	Production/activation rate of astrocytes by Aβo	1.793/day estimated
λ_A_O	Production rate of A βO	5×10−2/day estimated
λ H	Production rate of HMGB-1	3×10−5/day estimated
λ MF	Production/activation rate of microglias by NFT	2×10−2/day estimated
λ MA	Production/activation rate of microglias by astrocytes	2.3×10−3/day estimated
λ_M1T_β	Rate of M 1→M 2	6×10−3/day estimated
λ_M^_1T_β	Rate of M^_1→M^_2	6×10−4/day estimated
λ_T_βM	Production rate of TGF- β by M	1.5×10−2 day −1 [56, 99]
λ_T_βM^	Production rate of TGF- β by M^	1.5×10−2 day −1 [56, 99]
λ_T_αM1	Production rate of TNF- α by M 1	3×10−2 day −1 estimated
λ_T_αM^_1	Production rate of TNF- α by M^_1	3×10−2 day −1 estimated
λ_I_10M_2	Production rate of IL-10 by M 2	6.67×10−3 day −1 [47, 90]
λ_I_10M^_2	Production rate of IL-10 by M^_2	6.67×10−3 day −1 [47, 90]
λ PA	Production rate of MCP-1 by astrocytes	6.6×10−8 day −1 estimated
λ_PM_2	Production rate of MCP-1 by M 2	1.32×10−7 day −1 estimated
θ	M 2/M 1 effectivity in clearance of Aβo	0.9 estimated
α	Flux rate of macrophages	5 estimated
β	Proinflammatory/anti-inflammatory ratio	10 estimated
γ	I 10 inhibition ratio	1 estimated

Table 3

Parameters’ description and value

Parameter	Description	Value
d_Aβi	Degradation rate of Aβi	9.51/day [82]
d_Aβo	Degradation rate of Aβo	9.51/day [82]
d_AβoM	Clearance rate of Aβo by microglia	2×10−3/day estimated
d_AβoM^	Clearance rate of Aβo by macrophages	10−2/day estimated
d τ	Degradation rate of tau proteins	0.277/day [88]
d_F_i	Degradation rate of intracellular NFT	2.77×10−3/day estimated
d_F_o	Degradation rate of extracellular NFT	2.77×10−4/day estimated
d N	Death rate of neurons	1.9×10−4/day estimated
d NF	Death rate of neurons by NFTs	3.4×10−4/day estimated
d NT	Death rate of neurons by TNF- α	1.7×10−4/day estimated
d_N_dM	Clearance rate of dead neurons by M	0.06/day estimated
d_N_dM^	Clearance rate of dead neurons by M^	0.02/day estimated
d A	Death rate of astrocytes	1.2×10−3 day −1 estimated
d_M_1	Death rate of M 1 microglias	0.015 day −1 [47, 74]
d_M_2	Death rate of M 2 microglias	0.015 day −1 [47, 74]
d_M^_1	Death rate of M 1 macrophages	0.015 day −1 [47, 74]
d_M^_2	Death rate of M 2 macrophages	0.015 day −1 [47, 74]
D_A_O	Degradation rate of A βO	0.951/day estimated
d H	Degradation rate of HMGB-1	58.71/day [95]
D_T_α	Degradation rate of TNF- α	55.45 day −1 [47, 74]
d_T_β	Degradation rate of TGF- β	3.33×102 day −1 [56, 99]
d_I_10	Degradation rate of IL-10	16.64 day −1 [47]
d P	Degradation rate of MCP-1	1.73 day −1[47, 74]
R 0	Initial inflammation by ROS	6 estimated
M 0	Monocytes concentration in blood	5×10−2 estimated
N 0	Reference density of neuron	0.14 g/c m 3 estimated
MG0	Source of microglia	0.047 g/c m 3 estimated
A 0	Reference density of astrocytes	0.14 g/c m 3 estimated
K¯_Aβo	Michaelis-Mention coefficient for Aβo	7×10−3 g/ c m 3 estimated
K¯_N_d	Michaelis-Mention coefficient for N d	10−3 g/ml estimated
K_I_10	Half-saturation of IL-10	2.5×10−6 g/ c m 3 estimated
K_T_β	Half-saturation of TGF- β	2.5×10−7 g/ml [90]
K M	Half-saturation of microglias	0.047 g/ml estimated
K_M^	Half-saturation of macrophages	0.047 g/ml estimated
K_M_1	Half-saturation of M 1 microglias	0.03 g/ml estimated
K_M_2	Half-saturation of M 2 microglias	0.017 g/ml estimated
K_M^_1	Half-saturation of M^_1 macrophages	0.04 g/ml estimated
K_M^_2	Half-saturation of M^_2 macrophages	0.007 g/ml estimated
K_F_i	Half-saturation of intracellular NFTs	3.36×10−10 g/ml [89]
K_F_o	Average of extracellular NFTs	2.58×10−11 g/ml estimated
K_A_O	Average of of A βO	1×10−7 g/ml estimated
K P	Half-saturation of MCP-1	6×10−9 g/ml estimated
K_T_α	Half-saturation of TNF- α	4×10−5 g/ml estimated

We next observe that neurons are dying at approximately the rate of 5% a year, which was one of our important assumptions that was based on clinical data. We also note that, as the disease progresses, the plaque of A β peptides, , and the soluble A β oligomers, A , are increasing; reaches the level of 7×10−6 g/ml, in agreement with clinical data [57], and the assumed average of A concentration, , is indeed in good approximation to the average of the profile of A in Fig. 2. The assumed average of the F concentration, , is also in good agreement with the average of the profile of F in Fig. 2.

We note that N nearly stabilizes over time, at the level assumed in “Appendix," which means that, over time, macrophages and microglias clear debris of dead cells at nearly the same rate at which neurons are dying. Hence becomes very small over time, resulting in significant decline in extracellular NFT, while intracellular NFTs (F ) maintain a comparatively high level.

We finally note that the density of activated astrocytes is slightly increasing in agreement with a mouse model [58] which reports that astrocytes become increasingly prominent with the progression of the disease. The increase in A causes P also to increase, and the average of P is approximately equal to our estimate of K in S.I.

Anti-Alzheimer drugs

Until now, all clinical trials aimed to develop drugs that can cure AD have failed. There are currently no drugs that can prevent, stop or even delay the progression of Alzheimer’s disease, and there are many ongoing clinical trials. According to the 2016 Alzheimer’s Disease Facts and Figures, and the National Institute of Aging, if no cure is found, by 2050 the number of alzheimer’s patients in the U.S. will reach 15 millions and the cost of caring for them will exceed $ 1 trillion annually.

Avenues for AD therapies include prevention of build up of plaque (anti-amyloid drugs), preventing tau aggregation, and reducing inflammation. Clinical trials are concerned with both safety and efficacy. Here we shall use our mathematical model to conduct in silico trials with several drugs, addressing only the question of efficacy.

Treatment for AD causes changes in the densities of cells and concentrations of cytokines. In order to determine the efficacy of a drug, we should observe (i) to what extend it decreases the death rate of N, since slowing the death of neurons will improve cognition of patients; and (ii) to what extend it decreases , since A β aggregation mediates rapid dysfunction of synaptic plasticity and dendritic channels thereby causing memory loss [36-39].

TNF- inhibitor

Since TNF- α is implicated in generating neurotoxicity which leads to death of neurons, TNF- α inhibitor (etanercept) has been considered as a drug for Alzheimer’s patients [59]. In 2015 clinical trials phase 2 [60] the drug has shown some favorable trends but with “no statistically significant changes in cognition.” Since there were no serious adverse events, it was suggested that a larger, broader group needs to be tested before recommending etanercept for use for general Alzheimer patients.

We shall apply our model to determine how this TNF- α inhibitor affects AD patients. We use the following procedure:

Run the model for 300 days in order to ensure that AD has been diagnosed in patients;

Apply continuous treatment by the drug from day 300 until the end of 10 years.

During treatment, the effect of the drug is to replace Eq. (17) for TNF- α by the equation

where f is proportional to the amount of etanercept. We note that since etanercept is a soluble TNF receptor fusion protein, it stabilizes TNF- α [61] and thus TNF- α is diminished at rate f T . The red profiles in Fig. 3 show the result of the treatment with , compared to no treatment.

Fig. 3

Anti-TNF- α drug (red), etanercept, with ; TGF- β injection. (light-blue) with . Dark-blue color corresponds to no treatment, and where several profiles nearly coincide, they are all colored by light-blue. All the other parameters are as in Tables 2 and 3

TGF- injection

TGF- β is an anti-inflammatory cytokine which induces phenotype change from proinflammatory to anti-inflammatory macrophages. It was suggested that TGF- β mitigates AD pathology [29-34].

We note that the effect of T injection is to decrease M 1 and (see Eqs. (11), (13)), which results in a decrease in T (by Eq. (17)) and hence in a decrease in neuronal death rate. To model the treatment by injection of TGF- β we replace Eq. (15) for T by the equation

where g is proportional to the amount of injected TGF- β. In steady state, T maintains the level of , while its degradation rate is . Hence the source of T in steady state is . We take g to be 10 times this source, that is . We then follow the same treatment procedure for TNF- α inhibitor. The light-blue profiles in Fig. 3 show the results of the treatment, compared to no treatment.

Anti-A drugs

There are several drugs in Phase 3 clinical trials that aim to reduce the effect of A β aggregation [62]. Among them is aducanumab, which is thought to be microglia-mediated phagocytosis and clearance of A β [63, 64]. In our model, this drug will cause a decrease in the concentration of soluble , by replacing Eq. (2) by the equation

where h is proportional to the amount of the dozing level; we take h=10.

Figure 4 shows the efficacy of several drugs in terms of N and . The lowest curve in Fig. 4 a, and the highest curve in Fig. 4 b, correspond to the case where the curves of no treatment and several other drugs coincide; these drugs have negligible efficacy. Following the same treatment procedure as in the case of TNF- α inhibitor, Fig. 4 shows no efficacy of aducanumab in terms of N but significant efficacy in terms of in comparison to no treatment and to treatments by TNF- α inhibitor and TGF- β injection.

Fig. 4

Treatment with etanercept (decreasing T degradation rate by 10 fold), TGF- β injection (increasing its constitutive source by 10 fold), aducanumab (increasing the clearance rate of by 10 fold), and bindarit (increasing MCP-1 natural degradation by 10 fold). In a, the profiles of no treatment, bindarit and aducanumab coincide. In b, no treatment and bindarit coincide. The lowest curve in Fig. 4 a, and the highest curve in Fig. 4 b, correspond to the case where the curves of no treatment and several other drugs coincide; these drugs have negligible efficacy

MCP-1 inhibitor

Bindarit was shown to inhibit CCL2 (MCP-1) in brain tissue [65]. Hence it decreases (by Eq. (13)), which results in a decrease in T (by Eq. (17)) and thus also in a decrease in neuronal death rate. Bindarit was also reported to inhibit A β-induced neuronal death in vitro [66]. Hence it has a therapeutic potential in the treatment of neuroninflammatory/neurodegenerative diseases like AD [66]. We can conduct in silico trial with bindarit by revising Eq. (18) for P, replacing it with the equation

with k=10 Following the treatment procedure as in the case of of TNF- α inhibitor, Fig. 4 shows no efficacy of the drug in terms of N and in comparison to no treatment.

Methylthiomnium chloride (MTC) is the first identified tau aggregation inhibitor currently in Phase 3 trial [27]. In our model the drug will cause a decrease in the production of tau proteins and in their ability to turn into NFT. We model this by multiplying the production terms λ and λ by 1/10. Following the procedure as in case of TNF- α inhibitor, we found that the drug has almost negligible efficacy (not shown here).

Combination therapy

The results of Fig. 4 suggest that a combination therapy with etanercept (TNF- α inhibitor) and aducanumab (anti-A β drug), under the same ’10-fold’ amount, could both slow the death rate of neurons and decrease the growth of A β in a significant way. Figure 5 shows the dynamics of N and under such combination of drugs with different proportions of fold numbers: (etanercept,aducanumab)=(0,0) (no drugs), (10,5), (20,10), (30,15), (40,20) and (50,25). The reduction in the death of neurons, after 10 years, compared to the case of no drugs, is 3.8, 5.2, 6.4, 7.9 and 9.2%, respectively, and the respective reduction in the concentration of A β is 21, 32.2, 43.6, 53.9 and 64.1%.

Fig. 5

Combined treatment with etanercept fold number f and aducanumab fold number h for several values of (f, h)

We next consider combined therapy for any value of etanercept (f) and aducanumab (h). We define the N-efficacy of (f, h), E (f,h), to be where the density N is computed at the end of 10 years. Similarly we define the anti- efficacy by where the concentration is computed also at the end of 10 years.

Figure 6 is an efficacy map of the combined therapy with f in range of (0,50) and h in the range of (0,25). For any pair (f,h) the color columns in Fig. 6a and 6b show the efficacy for N and anti-.

Fig. 6

Efficacy maps. Etanercept (with fold number f) varies along the horizontal axis, and aducanumab (with fold number h) varies along the vertical axis. The column vector indicates the efficacy of treatment for any pair (f,h): a N-efficacy; b Anti- efficacy

We see that the efficacy of the combined therapy is very small if f<20 or h<10, and it increases sharply with f and h in the region where {40

From Fig. 6 we see that anti-A β antibody decreases the external concentration of A β () with efficacy less than 0.5 (h=20, f=0). Higher efficacy requires T inhibitor (h=20, f=20) which will protect neuron from death and prevent astrocytes activation, and thereby reduce . This result can be explained by our assumptions in Eq. (2) where we neglected the production of by live neurons and the increase of by ROS.

The PK/PD literature employs the concept of combination index (ϕ) in order to assess the level of synergy between two drugs [67]. This concept was used in simulations of several diseases (e.g. cancer and microbial diseases) in order to determine optimal dosage regimens [67-69]. Since in our AD model it is not clear how to define ϕ, and no data are available to evaluate ϕ, we shall, instead, introduce the following concept, for example in the case of etanercept and aducanumab:

We say that these two drugs at concentrations f and g have positive synergy with respect to N if E (f,g)>E (2f,0) and E (f,g)>E (0,2g). We accordingly define the synergy index σ =σ (f,g) by

Thus, σ >1 means positive synergy and σ <1 means negative synergy. The above definition depends on the doses f, g. If σ is large then the combination therapy at total amount f+g is much more effective than a single therapy, at the total same amount, in reducing the death rate of N. If σ <1 then a single drug is preferable. Similarly one can define the synergy index with respect to . Figure 7 shows the synergy index σ for (f, g) in the range 0aducanumab (h). Furthermore, given a total amount A of the combined drugs, so that f+h=A, the synergy increases as f/g increases. This suggests that in an optimal regimen f should be significantly larger than h, provided negative side-effects are discounted.

Fig. 7

Synergy map for combination therapy with etanercept (f) and aducanumab (h)

The synergy map for is similar to that of σ (not shown here), and so the synergy increases when f/g is increased.

From Fig. 5 we see that although the amyloid level are controlled, cell death levels do not decrease significantly. This may suggest that other combinations of drugs may target complimentary pathways more efficiently. For example, it was suggested in [70] that Amyloid β and tau combine to induce neuron into cell cycle, which leads to cell death; accordingly, one could explore using anti-A β and anti tau aggregation in combination therapy.

Sensitivity analysis

Sensitivity analysis on the model parameters can support the robustness of the simulation results. But it can also suggest what drugs do not work and what drugs are more likely to work. We conducted sensitivity analysis on parameters associated with production and removal rates of , death rates of N, and production rates of TNF- α, TGF- β and MCP-1: where we varied λ , λ , d , d between and twice their value in Tables 2 and 3, and varied γ, δ, ξ, ε between and 2.

Following the sensitivity analysis method described in [71], we performed Latin hypercube sampling and generated 2000 samples to calculate the partial rank correlation coefficients (PRCC) and p-values with respect to the density of N and with respect to the concentration of at time t=10 years. The results are shown in Fig. 8. All the p-values were less than 0.01. A positive PRCC (i.e. positive correlation) for N means that an increase in the parameter will increase the number of live neurons. A negative PRCC for N means that an increase in the parameter will decrease the number of live neurons. Similarly, a positive (negative) PRCC for means that an increase in the parameter will increase (decrease) the concentration of . Thus, for example, we see that d and d are negatively correlated to N and positively correlated to . This is not surprising since, with an increase in d and d , more neurons die (so N decreases) and as a result more A β emerge from the increasingly dying neurons, thus raising the concentration of . The fact that the correlation coefficients of d are significantly larger than the correlation coefficients of d , suggests that a drug which blocks TNF- α would be more effective than a drug which clears the F . The other PRCC values can also be seen to be consistent with the model dynamics.

Fig. 8

The PRCC values of parameter for sensitivity analysis

We observe that ε is negatively correlated to . Indeed, if ε is increased, more are cleared out (by Eq. (2)). To see how this affects N we note that if is decreased then A decreases (by Eq. (9)) and correspondingly M 1 decreases (by Eq. (11)), and then T decreases (by Eq. (17)); so we may expect N to increase, but perhaps not much, since we have ignored other indirect interactions from the model. From Fig. 8 we see that ε is indeed positively correlated to N but the correlation is small. The correlation levels of ε with respect to N and suggest that an anti-A β drug, like aducanumab, will have some benefits in reducing Amyloid β, but little benefit in reducing death of neurons. This is also seen from Fig. 4.

Conclusion

AD is an irreversible progressive neuroninflammatory/neurodegenerative disease that destroys memory and cognitive skills. Currently there is no drug that can cure, stop, or even slow the progression of the disease. Life expectancy at diagnosis is 10 years, and, at death, 50% of the brain neurons have already died. AD patients show abnormal aggregation of beta-amyloids () and neurofibrillary tangles (NFTs) of hyperphosphorylated tau proteins. NFTs destroy microtubles in neurons, which results in neurons death. Soluble oligomers activate microglias (the resident macrophages in the brain), thereby initiating inflammatory response. Additionally, peripheral macrophages, responding to cue from MCP-1 produced by astrocytes, are attracted to the brain and increase the inflammatory environment, which is harmful to neurons.

Figure 1 is a schematic network of AD: it includes neurons, astrocytes, microglias, peripheral macrophages, β-amyloids, tau proteins, and several cytokines involved in the cross-talk among the cells. In the present paper, we developed a mathematical model of AD based on Fig. 1. The model can be used to explore the efficacy of drugs that may slow the progress of the disease. We conducted several in silico trials with several drugs: etanercept (TNF- α inhibitor), injection of TGF- β, aducanumab (Anti-A β drug) and bindarit (MCP-1 inhibitor). We found that at ’10-fold’ level, etanercept has the largest efficacy in slowing death of neurons, while aducanumab has the largest efficacy in reducing the aggregation of , although these efficacies were quite small. Based on these findings we propose that clinical trials should use a combination therapy with etanercept (f) and aducanumab (h). In Fig. 6 we developed efficacy maps for any combination therapy with 0patient. When these limits become better known, one could then proceed to determine the optimal combination of etanercept and aducanumab for slowing the progression of AD.

The mathematical model developed in this paper depends on some assumptions regarding the mechanism of interactions involving amyloid, tau and neunofilaments in AD. There are currently not enough data to sort out competing assumptions. Hence the conclusion of the paper regarding combination therapy should be taken with caution.

Our mathematical model focused on the progression of AD in terms of neurons death and amyloid β aggregation. But dendritic pathologies also play an important role in the disease. Dendritic abnormalities in AD include dystrophic neuritides, reduction in dendritic complexity and loss in dendritic spines [36, 37]. In particular, A β plaques affect dendritic channels, and NFT mediates synaptic dysfunction [36-39]. Recent studies also begin to address white matter degeneracy that could help identify high risk of AD [72].

Appendix

Parameter estimation

In the sequel, in an expression of the form in the context of activation, the half-saturation parameter K is taken to be the steady state of the species X provided X tends to a steady state. Hence in a steady state equation this factor is equal to . If X does not tend to a steady state then the parameter K will be taken to be the estimated average of X over a period of 10 years, the average survival time of AD patients [73]. In an expression of the form (where γ=γ(X)) in the context of inhibition, K is again the half-saturation of X, so that in steady state the inhibition is 1/(1+γ). If cells Y phagocytose species X, then the clearing rate is proportional to where the Michaelis-Menten constant depends only on the ‘eating capacity’ of Y, so has no relation to the half-saturation of X.

Diffusion coefficients

The diffusion coefficient of proteins (Y) are proportional to , where M is the molecular weight [74]. Accordingly, we have the following relation [75]: where M and D are the molecular weight and the diffusion coefficient of VEGF. Since D =8.64×10−2 c m 2 day −1[76], M =24 kDa [76], M =8.9 kDa [77], kDa [78], kDa [79], kDa [46], and M =29 kDa [80], we get D =1.20×10−1 c m 2 day −1, c m 2 day −1, c m 2 day −1, c m 2 day −1 and D =8.11×10−2 c m 2 day −1.

Molecular weight of A β is 24 kDa [81], so in soluble state its diffusion coefficient would be 8.64×10−2 c m 2 day −1. We assume that soluble oligomer A β O has a smaller diffusion coefficient, namely, c m 2 day −1.

Eq. (1)

By [82], the half-life of is 1.5–2 h in mice. Hence =9.51 /day. Membrane proteins APP shed amyloid β, some end up inside the cell and some outside the cell. We assume that in healthy steady state , however the simulation results do not change appreciably if we take . According to [57], the density in brain-gray matter of is approximately 1000 ng/g in control and 7000 ng/g in AD. Hence, from the steady state of Eq. (1) in a healthy normal case, g/ml and = 9.51×10−6 g/ml/day. From the steady state of Eq. (1) in AD and Eq. (21) we then get that R 0=6.

The brain has 75% water and 60% of its dry matter is fat. We assume that the average density of brain tissue is 1 g/c m 3. The human brain has 100 billion neurons, and its weight is approximately 1400 g, so its volume is approximately 1400 ml. Hence its neurons number density is 7×107 neurons/ c m 3. The diameter of neurons is 16 μ m [83]. Accordingly, we estimate the volume of 1 neuron to be 2×10−9 c m 3, and the neurons density is then 7×107×2×10−9 g/c m 3, that is N 0=0.14 g/c m 3.

Eq. (2)

The number of neurons is three times the number of microglia [55], hence g/ml.

By [16] an astrocyte produces much less A β than a neuron, so we take .

Microglias are the first responders to NFTs and A βO. Peripheral macrophages arrive later, and their immune response may perhaps exceed that of microglia, but this is currently not known [12, 84]. We assume that in steady state the microglias density M and the peripheral macrophages density are equal, so that g/ml. Motivated by the inflammatory immune attack in AD [85], we assume that, in steady state, the proinflammatory macrophages exceed the anti-inflammatory macrophages, and that proinflammatory peripheral macrophages exceed the proinflammatory microglias. Thus, in steady state, , M 1>M 2 and , and we take , , , .

Activated microglias are poorly phagocytic for A β compared to peripheral macrophages [6]. Accordingly we take

Taking /day, we then have

We assume that and M 1 are more effective than and M 2 in clearing A β, and take θ=0.9.

We assume that survival time of patients with AD is 10 years, and that at the end-stage 50% of their neurons have died [73]. Hence, the death rate of N is /day.

By [57], g/ml. We assume that the clearance of by macrophages and microglias is nearly unlimited (i.e., it is almost linear in ) by taking g/ml. To estimate λ , we first consider the steady state of Eq. (2),

To estimate the average of , we use the equation so that

The values of for 500

The estimate of λ was based on the steady-state assumption in Eq. (2). However, in AD the A β peptides are continuously aggregating, so that the steady state assumption needs to be revised. We do this by increasing the value of λ : we take λ =2×4×10−9= 8×10−9 g/ml/day, and then λ =8×10−10 g/ml/day.

The number of astrocytes is approximately equal to the number of neurons [86, 87], hence A 0=N 0=0.14 g/ml.

Eq. (3)

Half-life of tau proteins is 60 hours [88]. Hence =0.277/day. Concentration of tau proteins is in healthy normal individuals is 137 pg/ml and, in AD, 490 pg/ml [89]. From the steady state of Eq. (3) in the healthy case, we have λ =d τ, where τ=137 pg/ml. Hence λ = 3.78×10−11 g/ml/day. Similarly, λ +λ R=d τ in AD, where τ=490 pg/ml. Hence we have λ R=8.1×10−11 g/ml, or λ =1.35×10−11/day.

Eqs. (4) and (5)

We assume that neurofibrillary tangles inside neurons are much more stable than tau proteins, taking /day. We also assume that extracellular NFTs do not degrade as fast as internalized NFTs, taking /day.

We also assume that 60% of the hyperphosphorglated tau proteins become neurofibrillary tangles. From the steady state of Eq. (4) we then have that . Hence λ =1.662×10−3/day.

Eq. (6)

It is not known whether the rate of death of neurons caused by NFT is larger or smaller than the death rate caused by T . We take d =2d , but the simulation of the model in the case where d =2d are very similar (not shown here). Assuming that at steady state of Eq. (6) the concentrations of F , T and I 10 are at half-saturation, we get , so that /day and /day. In particular, if γ=1 then d =2.4×10−4/day and d =1.7×10−4/day. We take g/ c m 3 (which is somewhat larger than the estimated half-saturation of I 10 in lung inflammation [47, 90]). We assume that in AD, 60% of hyperphosphorylated tau proteins (whose concentration in disease is 490 pg/ml [89]) are in NFT form, so that pg/ml= 2.94×10−10 g/ml. In [89] the concentration of tau protein was taken uniformly in the tissue of patients. We assume, however, that the concentration of NFT is higher inside neurons than outside neurons, and take g/ml, g/ml. From the steady state of Eq. (17) and the estimates of and (see under Eq. (17) below) we get T =4×10−5 g/ml, so that g/ml.

Eq. (7)

We take the half-life of astrocytes to be the same as the half-life of ganglionic glial cells, that is, 600 days [91]. Hence d =1.2×10−3/day. We assume that the activation of astrocytes is due more to TNF- α than to A β, and take . By the steady state of Eq. (7) we then get /day, and /day. Actually, in a mouse model of AD, the number of activated astrocytes is increasing [58]. So we compensate for this by increasing both and by a factor 1.1, taking /day and /day.

Eq. (8)

In mice experiments [92], macrophages phagocytosed apoptotic cells at rates that varied in the range 0.1–1.27/h. We assume that necrotic cells (and their debris) in human brain are phagocytosed by peripheral macrophages at rate /day. We also assume that microglia play a greater role in clearing necrotic neurons, and take =0.6/day. We also take g/ml.

Eq. (9)

We assume the degradation rate of A is much slower than that of , taking /day. The ratio of soluble A to total is approximately [93].

From the steady state of Eq. (9) we then get /day.

The estimate of was based on the steady-state assumption in Eq. (9). However, in AD the soluble A β oligomer is continuously increasing, following the increase in , so the steady-state assumption needs to be revised. We do this by increasing the above value of λ , taking the new value to be λ =5×10−2/day.

Eq. (10)

Concentration of HMGB-1 in neurons is 1.3 n g/m l [94], hence H=0.14×1.3 ng/ml= 1.8×10−10 g/ml. Half-life of HMGB-1 is 17 minutes [95], so that d =58.71/day. We assume that N stabilizes somewhere below 2.5×10−4 g/ml. From the steady state of Eq. (10), we then get λ =3×10−5/day.

Eqs. (11) and (12)

We take /day [47, 90]. Then, our assumption (under Eq. (2)) that suggests that β>1. We take β=10.

We take g/ml and α=5. In the absence of data, we take the production rate λ of macrophages by NFT to be the same as the production rate under stimulation by M. Tuberculosis in [90], namely, λ =2×10−2/day. We assume that production rate of macrophages by NFT is larger than the production rate by A , and take λ =2.3×10−3/day.

By [57] the concentration of A β in AD is 7×10−6 g/ml and, by [55], the ratio of A to is , so that g/ml.

We assume that more NFT reside within neurons than outside them, so that is smaller than . Recalling that g/ml, we take g/ml.

The coefficient is the rate by which TGF- β affects the change of phenotype from M 1 to M 2. In the case of infection in the lung by M. tuberculosis, under inflammatory conditions caused by the pathogen, /day [90]; we take it to be the same in the present case. We take g/ml, and g/ml.

Eqs. (13) and (14)

Peripheral macrophages immigrate into the brain of AD [96, 97]. We assume that, because of the BBB, the concentration of monocytes in the brain capillaries must be significantly higher than the concentration of peripheral macrophages already in the tissue. Recalling that in steady state g/ml, we take M 0=0.05 g/ml. The parameter α was estimated by 5, in order to make the asymptotic behavior of in the simulations agree with its assumed steady state of 0.047 g/ml (under Eq. (2)). When microglia cells are activated, they become either of M 1 or M 2 phenotype. But peripheral macrophages are initially biased toward phenotype rather than phenotype, since . We assume, in line with this bias toward , that the transition rate from into phenotype by TGF- β is at a smaller rate than the corresponding transition rate for microglias, that is, . We take /day.

Eq. (17)

Activated alveolar macrophages produce TNF- α at rate 4.86×10−3/day [47]. We assume that proinflammatory macrophages produce TNF- α at a larger rate (five fold), taking g/ml.

Eq. (18)

Astrocytes secrete MCP-1 [17-19] but activated anti-inflammatory microglias also secrete MCP-1. We assume that the production rate by astrocytes in larger than that by M 2, and take . MCP-1 concentration in initial stages of AD is 750 pg/ml [98]. Using the steady state equation

with P=6×10−9 g/ml and d =1.73/day [74], we get /day and λ =6×10−8/day [47].

Since A is increasing in time, also P is increasing in time. Hence the steady state assumption needs to be revised. We do it by increasing λ and by a factor 1.1, taking /day, and λ =6.6×10−8/day.