## Abstract

We present a model for chemotaxis, as applied to the aggregation of microglia in Alzheimer's disease. Using biological parameters found in the literature, testable thresholds are derived such that amyloid plaques are predicted not to form if conditions fall below a threshold. The model we use was developed by Luca *et al*. (Luca *et al*. 2003 *Bull. Math. Biol.* **65**, 693–730) and incorporates terms for both attractant and repellent signals. Our analysis can be applied to both two and three-dimensional spatial domains with application to any cell system involving chemotaxis.

## 1. Introduction

The formation of amyloid plaques is a key histopathological feature of Alzheimer's disease, a neurodegenerative disease (Selkoe 1997). The presence of β-amyloid causes neuronal cell death, but before this non-reversible event, other cell types contribute to the developing pathology. One of these cell types that is key to the pathology is microglia (Potter 1992). In response to β-amyloid there is a localized concentration of microglia (Itagaki *et al*. 1989) and an inflammatory response, which involves the local activation of microglia and the release of neurotoxins and factors that chemotactically attract more microglia from peripheral sites to the amyloid lesions (Rogers *et al*. 2002). Microglia also participate in the clearance of β-amyloid (reviewed in Rogers *et al*. 2002) and, therefore, help reduce local concentrations (Weldon *et al*. 1998). In an attempt to model the chemotactic response of microglial cells, Luca *et al*. (2003) have taken values from the literature for the density of microglia in normal and Alzheimer patient brain (Mackenzie *et al*. 1995), as well as concentrations and the likely decay of a probable attractants, such as IL-1β, and of a possible repellent, TNF-α, to give a prediction about amyloid plaque density and size. While biological data on Alzheimer's disease have been used to test the prediction, the Luca model and also our refinement presented here are best considered as a general model of chemotaxis.

We present here the further development of the model using the parameters presented by Luca *et al*. (2003). By applying nonlinear stability analysis to the model and parameters used by Luca *et al*. (2003) we derive a lower bound for a threshold concentration of microglial cells required to coincide with plaque formation.

The mathematical technique employed is somewhat akin to that of a generalized energy method in fluid dynamics, see e.g. Straughan (2004). The precise technique has much in common with the analyses for the Keller–Segel model of chemotaxis given by Horstmann (2001), Osaki *et al*. (2002) and Payne & Straughan (2004), although the variational analysis contained herein is new in a biological context as far as we are aware. We refer to Horstmann (2003, 2004) for detailed reviews of the Keller–Segel chemotaxis model.

## 2. Basic equations

Luca *et al*. (2003) study the chemoattraction–chemorepulsion model(2.1)In these equations *m*, *ϕ*, *ψ* represent, respectively, the cell density of microglia, the concentration of attractant (interleukin-1, IL-1β), and concentration of repellent (tumour necrosis factor-α, TNF-α). The summation convention applies to a repeated index in equation (2.1). In full generality *μ*, *Χ*_{1}, *Χ*_{2}, *D*_{1}, *D*_{2}, *a*_{1}, *a*_{2}, *b*_{1}, *b*_{2} may depend on *m*, *ϕ* and *ψ*. However, Luca *et al*. (2003) concentrate on the case where these functions are constant. They also restrict attention to one spatial dimension. In this work we likewise assume the above functions are constants although we allow the spatial domain to be two-dimensional. Our analysis can be applied in three-dimensions at the expense of an increase in technical detail. We briefly indicate in §4 how the extension to three-dimension proceeds.

Even though *Χ*_{1} and *Χ*_{2} are constant the chemoattraction and chemorepulsion terms are still nonlinear. Luca *et al*. (2003) developed a linear instability analysis and also investigated numerically simulated solutions. We here develop a fully nonlinear stability analysis. This leads to a nonlinear threshold below which aggregation (and, therefore, senile plaque formation) cannot occur.

Equation (2.1) is defined on a bounded domain and the boundary of *Ω* is denoted by *Γ*. Luca *et al*. (2003) introduce a length *L*, which is a typical dimension of a relevant domain in the brain. They do not give a value for this but instead when they non-dimensionalize the equations they use a spatial range *L*_{2}, which is the distance over which the repellent spreads during the characteristic time of decay. We employ exactly the same non-dimensionalization as Luca *et al*. (2003), details are given in eqns (6), (10), (11) of their paper. This leads to the following system of equations(2.2)The constants *A*_{1}, *A*_{2}, *ϵ*_{1}, *ϵ*_{2} and *a*^{2} are defined in Luca *et al*. (2003), but we include them below to make this article more readablewhere is the scale for the microglia density, and *L*_{1}, *L*_{2} are length-scales for attractant and repellent, respectively. The quantity is, in fact, the constant density of microglia in the steady state before any aggregation may occur.

The boundary conditions which hold are zero flux through *Γ*. This means that(2.3)where * n* is the unit outward normal to

*Γ*. Zero flux through

*Γ*for the conservation of microglia, equation (2.2)

_{1}signifiesbut due to equation (2.3) we must have(2.4)The constant steady state we are interested in is(2.5)although it has to be remembered that these are non-dimensionalized. In fact, must be consistent with its initial data,

*m*

_{0}(

*), so that where is the average value of the initial data*

**x***m*

_{0}over

*Ω*. (This point is discussed further in Payne & Straughan (2004) in the context of chemotaxis.)

## 3. Perturbation equations

The fully nonlinear perturbation equations from (2.2) are derived by setting , , . The perturbation quantity *m* has zero mean as a function of time, as may be seen by integrating equation (2.2)_{1} and using the boundary conditions (2.3) and (2.4), i.e. , where here with |*Ω*|=measure of *Ω*. The perturbation equations become(3.1)The functions *m*, *ϕ*, *ψ* satisfy the boundary conditions (2.3), (2.4), and equations (3.1) are to be solved subject to prescribed initial data(3.2)

## 4. Nonlinear solution decay

We now wish to develop a fully nonlinear stability analysis which will yield a meaningful threshold guaranteeing aggregation of microglia will not occur. Let ‖.‖ and (.,.) be the norm and inner product on *L*^{2}(*Ω*). We work with the solution measure(4.1)where *λ* and *λ*_{2} are positive constants we select optimally. The first term is the *L*^{2} measure of the perturbation to the microglia density, in the original variables this would be . Since , this is a natural measure to take and allows us to estimate how *m* decays to its steady state value . The *L*^{2} norms of ∇*ϕ* and ∇*ψ* are natural choices because equations (2.2)_{2,3} do not have right-hand sides which are gradients, like equation (2.2)_{1}. Decay of ∇*ϕ* and ∇*ψ* ensures the solution will not approach a state which is non-homogeneous in space. To proceed we multiply equation (3.1)_{1} by *m* and integrate over *Ω* using the boundary conditions (2.3) and (2.4) to derive(4.2)Next, multiply equation (3.1)_{2} by −Δ*ϕ* and (3.1)_{3} by −Δ*ψ* and integrate over *Ω* using the boundary conditions (2.3) to find(4.3)(4.4)We now add equations (4.2)+*λ*(4.3)+*A*_{2}(4.4), i.e. choose *λ*_{2}=*A*_{2} in equation (4.1). This has the effect of removing the *A*_{2} terms from the quadratic part of the analysis and means we shall not benefit from the stability effect of chemorepellent. However, we show in §5 that this is, in fact, the best we can achieve via the current route. (It is of interest to observe that an analogous feature is found in a nonlinear stability analysis for a solar pond, i.e. thermal convection in a layer of salty water which is heated and salted from below, cf. Straughan 2004, p. 238.)

With *A*_{2}=*λ*_{2} in equation (4.1) we find(4.5)We use the arithmetic–geometric mean inequality on the (∇*ϕ*, ∇*m*) term in equation (4.5) as follows:(4.6)where *α*>0 is a constant at our disposal. (This may seem a crude estimate to employ. However, we use a variational analysis in §5 to show it is, in an appropriate sense, optimal.)

Our analysis requires we keep a piece of ‖Δ*ϕ*‖^{2} in the dissipation term and so we write(4.7)for a constant 0<*ϵ*<1, and we use the inequality(4.8)on the *λϵ* part. The constant *λ*_{1}>0 and values for *λ*_{1}=*λ*_{1}(*Ω*) may be derived as indicated in Payne & Straughan (2004).

Combining equations (4.6) and (4.7) and the above use of inequality (4.8), in equation (4.5) we thus see that(4.9)where(4.10)In deriving equation (4.9) the cubic terms in equation (4.5) have been manipulated as follows:where we have integrated by parts and used the boundary conditions (2.3). A similar analysis holds for the *ψ* term.

To handle the cubic terms in equation (4.9) we use the Cauchy–Schwarz inequality to see that(4.11)and we then use the Sobolev inequality, valid for functions with mean value zero which satisfy equation (2.4) (given in e.g. Payne & Straughan 2004),where a value for *c*_{1} may be found in Payne & Straughan (2004). Upon estimating the cubic terms in this way we derive from equation (4.9)(4.12)where the dissipation term *D* is given by(4.13)From inequality (4.12), recalling the definitions of *E*(*t*) and *D*(*t*) we may show(4.14)where(4.15)Decay of *E*(*t*) follows from equation (4.14), provided we ensure two conditions hold. These are that the coefficients in *D* in equation (4.13) are positive and that *E*^{1/2}(0)<1/*K*. Details of how this leads to exponential decay of *E*(*t*) in time are not given here but may be found in e.g. Straughan (2004, p. 15). The conditions which guarantee *E*(*t*)→0 exponentially are of interest to us because this is what ensures aggregation of microglia cannot form. The analysis shows *m*→0 and ∇*ϕ*, ∇*ψ*→0 so the chemoattractant and repellent concentration must tend to a constant while the microglia density must return to its constant steady state value. We stress this analysis is fully nonlinear.

The conditional restriction *E*^{1/2}(0)<1/*K*, which restricts the size of the initial perturbations that can perturb the steady state in , , , is(4.16)The restrictions on the coefficients *r* and that of ‖∇*ϕ*‖^{2} in *D* require(4.17)The last condition implies we select *α* such that(4.18)Hence, we requireThe function *f*(*λ*) has a maximum at *λ*=(*ϵλ*_{1}+*a*^{2})/*a*^{4} and this leads to the restriction on *A*_{1},(4.19)Hence, we are guaranteed decay of *E*(*t*), and so non-aggregation of microglia, provided *A*_{1} satisfies inequality (4.19) and the initial perturbations satisfy restriction (4.16). The interpretation of these conditions is deferred until §6.

*Note*. If we perform an analogous analysis for , then the functional *E*(*t*) is not sufficient to yield decay. Instead one must employ a functional which includes terms in ‖Δ*ϕ*‖^{2} and ‖Δ*ψ*‖^{2} in addition to those present in *E*. The condition (4.19) is still achieved, although the conditional restriction (4.16) is more complicated.

## 5. Variational analysis

We now justify the choice of *λ*_{2}=*A*_{2} in §4 and the use of the arithmetic–geometric mean inequality in equations (4.6)–(4.9).

We commence with *E*(*t*) as defined in equation (4.1) with *λ*>0, *λ*_{2}>0 arbitrary and at our disposal. Instead of equation (4.5) we arrive at(5.1)where(5.2)(5.3)with *N* defined as in equation (4.11).

We proceed with a variational analysis which will produce a sharp result. Details of such variational calculations in a fluid mechanics setting may be found in the research monograph by Straughan (2004). The procedure is to define(5.4)where *H* is the space of admissible solutions. Then, from equation (5.1) deriveFor *Λ*≥1 we showand decay of *E* may be established. The decay criteria are now that *E*^{1/2}(0)<1/*K* and *Λ*≥1.

We solve the variational problem (5.4). The Euler–Lagrange equations are found to be(5.5)Suppose now the functions *m*, *ϕ* and *ψ* satisfy the equation Δ*ϕ*=−*q*^{2}*ϕ* where *q* is a wavenumber. A pattern formation function like this is discussed in a fluid mechanical context in e.g. Straughan (2004, p. 54); for example, if *Ω* is a rectangle then for a constant amplitude is a solution. Provided *q*≠0, we can then obtain from equations (5.5)(5.6)where *β*_{1}=(*A*_{1}+*λa*^{2})/2, *β*_{2}=(*λ*_{2}−*A*_{2})/2.

After some rearrangement equation (5.6) can be written as(5.7)For nonlinear stability we require *Λ*≥1, and so at the optimal value equation (5.7) requires(5.8)The parameters *λ*, *λ*_{2} are independent and are selected to optimize equation (5.8). The best value of *λ*_{2} clearly makes the second term on the left of equation (5.8) as small as possible and so *λ*_{2}=*A*_{2}. With this value equation (5.8) reduces to(5.9)Since Δ*ϕ*+*q*^{2}*ϕ*=0 condition (5.9) is the same as condition (4.18) when *ϵ*=1. The number *λ*_{1}≡*q*^{2}. Thus, we have shown by a variational approach that condition (4.19) is optimal, at least with the energy-like argument employed here.

## 6. Interpretation

While this section is highly speculative, especially since we do not know all the values of coefficients and we cannot be sure the model does have any bearing on Alzheimer's disease, we include it because we wish to show that our analysis may be useful in the context of the parameter values which are currently available.

The conditions to prevent aggregation into plaques of microglia are (4.16) and (4.19). Condition (4.16) depends on our choice of *λ* and *ϵ*. With the optimal value of *λ*=(*ϵλ*_{1}+*a*^{2})/*a*^{4} (as calculated before equation (4.19)), inequality (4.16) becomes(6.1)if we select *α*=(1+*ϵλ*_{1}/*a*^{2})^{−1} in accordance with equation (4.18).

Condition (4.19) is(6.2)(consistent with equation (6.1)). Given that *a*=*L*_{2}/*L*_{1} and in one dimension *λ*_{1}=*π*^{2}(*L*_{2}/*L*)^{2} this restriction on *A*_{1} meansThe number , and so we may reinterpret this inequality as a threshold on the microglia density below which aggregation cannot form. Note that now refers to the steady state microglia value before non-dimensionalization, i.e. the value appropriate to equations (2.1), not the value we use in equation (2.5)_{1} which is appropriate to equations (2.2). In this way we find the threshold(6.3)Using estimates from Luca *et al*. (2003) we findand we estimate a typical brain dimension *L*≈0.1 mm. Thus, if anything *L*_{1}≥*L*, which helps condition (6.3). The IL-1β decay rate calculated by Luca *et al*. (2003) is *b*_{1}=3×10^{−3} to 3×10^{−2} min^{−1}, while the reported chemoattraction coefficient *Χ*_{1} has value between 6 and 780 μm^{2} nM^{−1} min^{−1}. Using these values and others reported in Luca *et al*. (2003), we find a threshold value of from approximately 20.3 to 2.6×10^{4} nM pg^{−1}. To convert this to cells μm^{−3}, we assume a molecular weight appropriate to glia to be that of IL-1β (secreted by the glia), i.e. 17 kDa as given by Luca *et al*. (2003, p. 34; we realize this may be a crude estimate). With this figure we calculate the threshold for microglia density to be in the range(6.4)If we employ the molecular weight of glia to be that appropriate to β-amyloid plaques, i.e. 900 kDa, as suggested by Huang *et al*. (2000), then the threshold range of equation (6.4) increases to(6.5)This should be compared with a figure of microglia density in the brain of aor of anas calculated from the values given by Luca *et al*. (2003, p. 3).

At present equations (6.4) and (6.5) are not too useful because of the uncertainty of available values for the necessary coefficients. If is of the order of the lower value, , then this is of no value since it is smaller than that found in a healthy brain. If, however, the upper value of or 22.1×10^{−3} cells μm^{−3} is likely, then this could be significant in indicating a threshold below which aggregation to form plaques will not occur.

Of course, thresholds (6.4) and (6.5) have to hold simultaneously with the initial data restriction (6.1). The data restriction is not believed to be as significant as equations (6.4) and (6.5), however. In any case, at present no value is available for *Χ*_{2}, the chemorepellent coefficient, and so it is difficult to find a reliable value for the right-hand side of equation (6.1).

In future work, we intend analysing other factors such as the clearance rate of β-amyloid by microglia in the inhibitory response, cf. Ard *et al*. (1996), and other reaction constants, e.g. Hasegawa *et al*. (2002), in a modification of the Luca *et al*. (2003) model.

## Acknowledgments

We wish to thank two anonymous referees for critiscism which we have incorporated. This work is supported by NIH grant (NS42803) to R.A.Q.

## Footnotes

- Received November 22, 2004.
- Accepted March 23, 2005.

- © 2005 The Royal Society