The Chemicals

      Chemicals currently represented in the simulation include Beta-amyloid protein, interleukins IL-1B, IL-6 and tumor necrosis factor, TNF. (The former is responsible for initiating an inflammatory response, the latter three are cytokines secreted by glial cells.) Chemicals are produced and taken up by cells, diffuse throughout the domain and influence cell behaviour as described above in The cells.

Diffusion

      The diffusion coefficient of a substance is computed from its molecular weight, (Goodhill, 1998, 1997) and corrected for the overall effect of tortuosity of brain tissue (Sykova, 1997; Nicholson and Sykova 1998), and for possible local effects that affect diffusion (described below). For S(x,y,t) the concentration of some chemical at point (x,y) at time, t, and given its diffusion coefficient at point (x,y), i.e. D(x,y), the following partial differential equation (PDE) describes diffusion of the chemical,

.
By solving this equation numerically, we determine the concentration of the chemical in grid space (i,j), denoted by S1i,j, based on the concentrations of the chemical in the previous time step, denoted by S0i,j.

      Both explicit and implicit finite difference methods for numerically solving the chemical transport PDE were tested in our simulations. The constraints of running a Java-based simulation online have meant that we are currently using an explicit method in order to save on memory.

Explicit Method

      In explicit methods, spatial derivatives use data from the preceding time step. The discretized PDE then has the form

.
It is straight forward to solve for S1i,j so that
.
Letting i range from 1 to m and j range from 1 to n, one needs to use the above formula mn times to find the resulting concentrations throughout the environment. Even if the diffusion coefficients, Di,j, were to change, locally, due to changes in the evolving state of the tissue, the method works with little change in the cost of computing S1i,j.

      A draw-back is that the above method is numerically stable only if the maximum diffusivity, D, is less than DX*DY/(4*dt) (for the case where DX = DY). Our current formulation sets DX = DY = 10 and dt = 0.0125 so that the maximum diffusivity is D = 2000 microns2 per minute. This stability constraint was used in selecting an appropriate micro time scale.

Implicit Method

      The discretized form of the PDE in an implicit method typically has the form:

.
The quantity S1i,j in the above equation in expressed in terms of the unknown quantities S1i-1,j, S1i+1,j, S1i,j-1 and S1i,j+1 to determine. As i ranges from 1 to m and j ranges from 1 to n, the above system forms mn linear equations with mn unknowns. Finding the value of all the S1i,j requires solving such a system or, equivalently, computing the inverse of the mn x mn matrix.

      The implicit method described above is the backwards Euler method. In the program, we use a generalized Crank-Nicholson scheme which averages the spatial derivatives between the old data and the new data. (The fact that one can weight the average is what makes it generalized.) In our current formulation, we invert an mn x mn matrix. This requires a huge computational overhead, but as long as the diffusion coefficients do not change, this only needs to be done once. If the diffusion coefficients change (as is currently the case), this method should be avoided. However, in such a case, another possibility is to use an Alternating Direction Implicit (ADI) method: this uses the fact that diffusion in each direction leads to a tridiagonal matrix which can be solved on its own, without building and inverting the entire mn x mn matrix. (There are also shortcuts for inverting a tridiagonal matrix which reduce the computational overhead.)

      With an implicit method, one rarely has to worry about stability. Thus, updates can take place on a macro time scale. In practice, the consideration of numerical accuracy generally restricts one to smaller time steps nevertheless. Further, the savings due to fewer updates relative to the explicit schemes may be outweighed by the additional computation power needed to invert matrices. One must have the memory to store the matrices, and when multiple inversions are necessary, they carry a cost of (mn)3 which can be enormous. It is for these reasons that we currently use an explicit scheme.

Secretion

      Currently, the production of chemicals by cells is assumed to occur at a constant rate, r. Thus, if a cell secretes some chemical, S, into a grid space, the rate of change of S in that grid space is given by the equation

.
The chemicals diffuse in the domain, and are absorbed by other cells.



Amyloid Protein

      Amyloid protein occurs in soluble and fibrous forms. In the soluble form, it is secreted by (live) neuronal tissue, diffuses in the tissue, interacts with amyloid fibers, and is removed by microglia.

      A source of soluble amyloid secreted from neurons in the center of the environment forms the initial "stimulus" at the begining of a simulation. Amyloid fibers are placed randomly in the domain with frequency based on the initial fiber occupancy parameter, p. (Initially, for each grid space, a random variable uniformly distributed between 0 and 1 is generated; if this value is less than p, a random concentration of fibers uniformly distributed between 0 and MAXFIBERS is placed in the grid space.)

Secretion

      When a neuron's internal concentration of IL-1B exceeds the source triggering level the Neuron secretes soluble amyloid protein at its grid location. The rate of secretion, r, is determined as a function of the source concentration parameter, I. Specifically,

r = I (s/c) d
where s is the concentration of IL-1B within the neuron, c is the maximum IL-1B absorbed parameter, and d is a programmer defined constant (currently set to 10) that was used to make older runs comparable to the new ones. Secretion occurs on the micro time scale so that the addition of soluble amyloid in the grid space is
DS = dt*r.
However, the maximal allowable soluble amyloid concentration is given by the source concentration parameter, and its level is capped at this ceiling.

Diffusion

      Soluble amyloid diffuses according to the diffusion properties discussed above. Initially, the diffusion coefficient equals the diffusivity parameter throughout the region. However, in regions where astrocytes are activatedm astrocyte blocking may partially seal off certain regions and thus reduce the diffusion coefficient of chemicals in the given region.

Transition of Soluble Amyloid to Fibers

      Exchange between soluble and fibrous amyloid occurs in two ways. Fibers already present in any site can grow in the presence of soluble amyloid (sol to fiber transition rate parameter governs this rate). Fibers in the surrouding grid spaces can elongate into neighboring sites. In either case, the concentration of soluble amyloid in the grid space must exceed the critical sol-AB for fibers parameter, for fiber growth to occur.

      When fibers are present, growth occurs at a rate proportional to the product of the average surrounding fiber concentration, F (a weighted average with the center fibers counting more than the immediately surrounding fibers by a programmer defined constant, WEIGHT -- currently, WEIGHT=2.0, meaning that fibers in the center count twice as much as surrounding fibers), and the difference between the concentration of soluble amyloid, s, within the grid space and the critical sol-AB for fibers parameter, b. If f is the concentration of fibers in the grid space, then

where R is the sol to fiber transition rate parameter. Because fiber transitions happen on the macro time scale, the change in fiber concentration, Df, can be calculated as
Df = DT*R*(s-b)*F.
There are two restrictions on the size of Df. It can be no larger than s, otherwise one would be trying to change more soluble amyloid to fiber than is present. As well, Df + f cannot exceed MAXFIBERS. If either case occurs, Df is reduced to be the maximum concentration that still passes both criteria.

      In the case where no fibers are present in the grid space under consideration, the grid space may still gain fibers if no astrocytes occupy the space and the concentration of soluble amyloid exceeds the critical sol-AB for fibers parameter. This is done through the process of fiber nucleation. Fiber nucleation can depend on (1) surrounding grid spaces having fibers or (2) purely on the amount of soluble amyloid.

      In the first case, fiber nucleation depends on the weighted average of concentration of fibers surrounding the grid space, F, the maximum fiber concentration, MAXFIBERS, and the new fiber nucleation effectiveness parameter, n. These terms are combined to give a "probability" that the site undergoes nucleation, p = n*(F/MAXFIBERS). Every macro time step, nucleation is tested for via a Monte Carlo technique. Specifically, a uniformly distributed random number between 0 and 1 is generated. If its value is less than p, then nucleation occurs and the rate of fiber growth is determined as described in the preceding paragraph.

      In the second case, fiber nucleation depends on the concentration of soluble amyloid. The new fiber nucleation effectiveness parameter, n, is divided by the maximum value that it can attain from its slider, to determine a "probability" of nucleation, p = n/max. Every macro time step, DT, nucleation is tested for via a Monte Carlo technique. Specifically, a uniformly distributed random number between 0 and 1 is generated. If its value is less than p, then nucleation occurs. In this case, fiber growth is

DF = DT * e * (S - M)2
where e is the sol to fiber transition rate, S is the concentration of soluble amyloid and M is the critical sol-AB for fibers. Currently, this case has a far better chance for fiber nucleation to occur than the first case (p is much greater in this case). It may be worthwhile to decrease p for this case so that it behaves more like the first. One method for doing this would be to allow p to behave like in the first case and simply set F to be some very small value like 1, so that p = n/MAXFIBERS.



IL-1B

      IL-1B is secreted by microglia and absorbed by astrocytes and neurons. It is free to diffuse throughout the environment.

Secretion

      Microglia secrete IL-1B into the same grid space at which they are located when their internal concentration of soluble amyloid exceeds the triggering concentration parameter. The appropriate checks are made so that the microglia concentration is taken into account. The rate of secretion is the product of the IL-1B secretion rate, u, and the microglia concentration, q, so that

r = q u.
Because secretion occurs on the micro time scale, the increase of IL-1B, DS, in the grid space of the microglial agent satisfies
DS = dt*q*u.
Currently, there are no limits to how much IL-1B can be in any one grid space.

Diffusion

      IL-1B diffuses according to the diffusion properties discussed above. Initially, the diffusion coefficient equals the diffusivity parameter throughout space. However, astrocyte blocking may reduce the diffusion coefficient in specific regions.



IL-6

      IL-6 is secreted by astrocytes and absorbed by neurons. It is free to diffuse throughout the environment.

Secretion

      Astrocytes secrete IL-6 into the same grid space at which they are located based on the amount of IL-1B they have in storage, s. Adjusting for the astrocyte concentration parameter, q, secretion occurs as long as s is greater than

q*kd*g,
where kd is the IL-1B receptor equilibrium constant for astrocytes and g is the triggering IL-1B equilibrium factor for IL-6. The rate of secretion is the product of the IL-6 secretion rate, u, and the astrocyte concentration, q, so that
r = q u.
Because secretion occurs on the micro time scale, the increase of IL-6, DS, in the grid space of the astrocyte agent satisfies
DS = dt*q*u.
Currently, there are no limits to how much IL-6 can be in any one grid space.

Diffusion

      IL-6 diffuses according to the diffusion properties discussed above. Initially, the diffusion coefficient equals the diffusivity parameter throughout space. However, astrocyte blocking may reduce the diffusion coefficient in specific regions.



TNF

      TNF is secreted by astrocytes and absorbed by neurons. It is free to diffuse throughout the environment.

Secretion

      Astrocytes secrete TNF into the same grid space at which they are located based on the amount of IL-1B they have in storage, s. Adjusting for the astrocyte concentration parameter, q, secretion occurs as long as s is greater than

q*kd*g,
where kd is the IL-1B receptor equilibrium constant for astrocytes and g is the triggering IL-1B equilibrium factor for TNF. The rate of secretion is the product of the TNF secretion rate, u, and the astrocyte concentration, q, so that
r = q u.
Because secretion occurs on the micro time scale, the increase of TNF, DS, in the grid space of the astrocyte agent satisfies
DS = dt*q*u.
Currently, there are no limits to how much IL-6 can be in any one grid space.

Diffusion

      TNF diffuses according to the diffusion properties discussed above. Initially, the diffusion coefficient equals the diffusivity parameter throughout space. However, astrocyte blocking may reduce the diffusion coefficient in specific regions.



previous: Cells