Cell Formulation

      The three types of cells in the simulation, astrocytes, microglia and neurons, all secrete and absorb chemicals. Astrocytes and microglia can move. It is possible to explore the effects of cell division and death of microglia, though these are omitted in the default setting.

Absorption

      In the current simulation, there are two types of absorption models, Michaelis-Menten kinetics and receptor kinetics.

Michaelis-Menten Kinetics

      Michaelis-Menten kinetics describe the rate of absorption as a function of the concentration of chemical present, S. Two parameters are necessary for determining the rate of absorption, the maximum absorption rate (call it k) and the concentration of chemicals at which the rate of absorption is half its maximum (call it h). Then the change in concentration S can be described by the following differential equation:

.

Receptor Kinetics

Case 1

      A slightly more detailed model for absorption kinetics tries to describe the process through receptor kinetics. Let r be the concentration of receptors for a given cell. Then by using the kinetics of this receptor, specifically the rate at which the chemical S binds to the receptor, kf, and the rate at which it unbinds, kb, one can formulate a differential equation through the laws of mass action to describe the binding kinetics. If one were to assume that half of the receptors are active (the proportion of active receptors actually varies), then the change in concentration S can be described by the following differential equation:

Because the receptor binding equilibrium constant, kd, is usually given rather than kf, it is useful to convert the equation into the following form:
.

Case 2

      However, the latest addition of the simulation relaxes the assumption that half of the receptors are bound. In this case, external chemicals, So, can be absorbed into cells to become internal, Si, by binding to receptors, r, as shown in the following reaction:

.
r in this reaction can be viewed like an enzyme so that this enzymatic reaction is fully described by the Michaelis-Menten kinetics shown above. By letting R be the total concentration of receptors, then we can use the equation above with
k = k2 R   and   h = (kb + k2)/kf
to describe chemical absorption by cells.


Microglia

      The initial microglia count dictates the number of microglial agents initially in the simulation. The microglia concentration designates how many microglia are represented per agent. Thus, if the intial count is 40 and the concentration is 10, then 400 microglia are represented in the simulation. By using the concentration parameter one can change the number of microglia in the simulation without changing the computation power. The microglial agent is placed (or centered) within a grid space. The number of agents centered in any one grid space cannot exceed the maximum density parameter.

      At initialization, a random pairing of (x,y) is generated. If there is room (meaning that the maximum density parameter will not be exceeded) to place a microglial agent in this grid space, a microlgial agent will be placed in the grid space. If there is no room, a new random pairing of (x,y) is generated and tested. This process continues until a number of microglial agents equal to the initial microglia count are in place.

Absorption of Amyloid Protein

      Microglia absorb soluble amyloid according to the Michaelis-Menten kinetics discussed above. The maximum sAB microglia uptake rate is k and the half max amyloid binding conc is h. Because this occurs on the macro time scale, the amount of soluble amyloid changed in the same grid space as the microglial agent, call it DS, is

DS = q*DT*k*S/(h+S)
where q is the microglia concentration and S is the concentration of soluble amyloid in the same grid space as the microglial agent. DS may be modified if the cell is nearing its capacity for soluble amyloid defined by the fatal sol-AB dosage parameter. If DS is great enough so that if it were added to the stored soluble amyloid in the microglium, the storage concentration would exceed the product of the capacity and density, then DS is decreased to the maximum concentration that can be stored without exceeding capacity. Having determined DS, the concentration of soluble amyloid in the same grid space as the microglial agent is decreased by DS and the concentration of soluble amyloid within the agent is increased by DS.

      Microglia can also ingest amyloid fibers at a constant rate. The programmer defines a fiber ingestion fraction, z (currently, z=0.1), that is multiplied with the maximum sAB microglia uptake rate, k to determine the constant rate of ingestion. Because this occurs on the macro time scale, the amount of fibrous amyloid, F, changed in the same grid space as the microglial agent, call it DF, is

DF = q*DT*z*k*F
where q is the microglia concentration. If DF is greater than F, then it is adjusted to equal F. There is no maximum to the amount of fiber a microglium can ingest.

Neutralization of Stored Soluble Amyloid Protein

      In addition to absorbing soluble amyloid, it is assumed that microglia can neutralize the concentration of soluble amyloid that it has stored. It does this at a constant rate defined by the sol-AB neutralization rate parameter. In this case, the change in concentration within the cell is simply the product of the microglia concentration, the neutralization rate and the macro time step. Of course, care is taken so that the concentration of soluble amyloid in storage is never negative.

Secretion of IL-1B

      Microglia secrete IL-1B. For details on how this is done, please see IL-1B under Chemicals.

Movement

      The motility time delay parameter determines how many macro time steps a microglial agent must wait before moving. If there is no delay, then the speed of movement is simply the length of a grid space, DX, divided by the macro time step, DT. In general, if m is the number of steps between movement, then the corresponding speed is v = DX/(DT+1).

      Each microglial agent has a counter which gets decremented every time step. When the counter reaches zero the cell attempts to move according to the rules specified below. Also at this time, the counter gets reset to the motility time delay value. When microglial cells are initialized, the counter takes any integer value between zero and the time delay, uniformly distributed. This allows for movement to be asynchronous, meaning that not all microglia will move at the exact same time.

      Movement of microglia can be hampered by the presence of amyloid fibers. If the concentration of fibers in the grid space of the microglial cell is F and the concentration of fibers at which a cell has a 50-50 chance of being stuck is called G (this is the half-sticking fiber level parameter), then the the probability, P, that a cell gets stuck is

P = F/(F + G).
This is tested in a Monte Carlo fashion meaning that a uniformly distributed random variable between 0 and 1 is generated and if its value is less than P, then the cell is stuck contingent on the concentration of soluble amyloid present. If the concentration of soluble amyloid, S, is great enough, then normally stuck cells are free to move depending on the concentration of soluble amyloid at which a cell has a 50-50 chance of being unstuck (call this H which is the half-sticking sAB level parameter). The probability, P, of a cell to remain stuck is
.
This is tested in the same Monte Carlo fashion.

      Microglia move chemotactically up the soluble amyloid gradient. Thus, the first step in cell movement is to determine which of the surrounding grid spaces has the greatest concentration of soluble amyloid. The grid spaces sampled can be labeled as

	0 1 2
	3 4 5
	6 7 8
where 4 is the current grid space of the microglial agent. In the case where the cell is on a boundary, periodic boundary conditions are assumed. After sampling the immediate neighborhood, the direction of greatest concentration is noted. The cell will move in this direction in a Monte Carlo fashion with the probability of winning given by the chemotactic sensitivity parameter. (Simply put, the cell will move in the direction of highest concentration with a probability equal to the chemotactic sensitivity parameter.) If the designated direction loses or all concentrations are the same so that there are no designated directions, then a direction 0-8 is chosen at random, uniformly distributed. Having finally decided on a direction of movement, the cell must check to see if room is available for movement. Room exists if no astrocytes are in the grid space and the addition of another microglial agent will not make the number of microglial agents in the grid space exceed the maximum density parameter. If there is no room, then the microglial agent will stay in its current grid space.

Cell Death

      If cell death is turned on, then one of two factors will determine whether a microglial agent will die. If the concentration of absorbed soluble amyloid ever reaches the fatal sol-AB dosage parameter, then the cell will die. A cell could also die based on old age. We assume that the probability that a cell dies by age T is based on a sigmoidal distribution,

where A is the typical cell duration. Solving the above equation for T allows us to determine the time of death based on the cumulative distribution function
.
By generating random values of P between 0 and 1, uniformly distributed, we can determine a corresponding random age of death. Thus, when the cell is created (or initialized), the time of death is known and this value is simply placed in a death counter (as T/DT) that gets decremented every macro time step. If the fatal soluble amyloid dosage is never reached before the death counter reaches zero, then death takes place when the counter reaches zero.

Mitosis (Cell Division)

      If cell mitosis is turned on, then four criteria determine whether mitosis can occur. These factors are the length of cell cycle, half-sticking fiber level, poison effects on mitosis, and chance of mitosis. All four criteria must be satisfied in order for mitosis to occur.

      The length of cell cycle determines the time at which mitosis can occur. Every microglial agent has a counter which has the cell's current age (in number of macro time steps since the cell was initialized). The length of cell cycle, call it C, is converted into number of macro time steps, M = C/DT. Whenever, the age counter is a multiple of M, mitosis can occur. Of course, this process can be changed, so that all microglial agents do not divide at the same time as follows: The length of the cell cycle can determine the mean of a normal distribution with the standard deviation another parameter to be added to the simulation. The timing for mitosis can then be handled much like cell death above, where each time a microglial agent is initiated, a time for cell mitosis can be randomly generated based on the normal distribution (or any other distribution). Cells can then keep track of a mitosis counter which gets decremented every macro time step. When the counter reaches zero, the other criteria are checked and the counter gets reset (stochastically based on the mitosis distribution like it was at initiation).

      Because amyloid fibers are deemed hazardous to microglia, presence of fibers may prevent mitosis. In our current formulation, if the concentration of amyloid fibers in the same grid space as the microglial agent is greater than the half-sticking fiber level, then mitosis is prevented. If the concentration of amyloid fibers is less than the half-sticking fiber level, then the other criteria are checked.

      The amount of soluble amyloid within the microglia (having been absorbed) also affects mitosis. The poison effects on mitosis parameter accounts for this. If the parameter equals one, then the cell is deemed poisoned if any soluble amyloid is stored inside the cell. If the parameter equals zero, then soluble amyloid does not act as a poison which prohibits mitosis. Values in between dictate a percentage of the capacity (the fatal sol-AB dosage parameter) at which the cell is poisoned. In this case, let E be the value of the poison effects parameter, S be the concentration of soluble amyloid within the cell, q be the microglia concentration and C be the capacity of soluble amyloid within a cell. The microglial agent is poisoned if S is greater than q(1-E)C in which case mitosis cannot occur. Otherwise, the other criteria are checked.

      In our current formulation, all the other criteria are deterministic (the condtions are either met or not). To add a stochastic component to the model, the chance of mitosis parameter is checked in a Monte Carlo fashion. Thus, if all the other criteria have been met, a uniformly distributed random variable between 0 and 1 is generated. If this number is less than the chance of mitosis parameter, mitosis occurs. Otherwise, nothing happens until the next time mitosis can be checked.

      If all the conditions for mitosis are met, then the program will attempt to add another microglial agent to the simulation. It first checks to make sure that there is enough memory (currently, the simulation can hold 1000 microglial agents). A new microglial agent is then initiated in the same grid space as the old one. The concentration of soluble amyloid in storage in the old agent is divided so that the old agent holds half and the new agent begins with half. If the addition of the new microglial agent would make the number of microglial agents exceed the maximum density parameter, then the new agent simply replaces the old one.



Astrocytes

      The initial astrocyte count designates the number of astrocyte agents initially in the simulation. The astrocyte concentration designates how many astrocyte are represented per agent. Thus, if the intial count is 100 and the concentration is 10, then 1000 astrocytes are represented in the simulation. By using the concentration parameter, one can change the number of astrocytes in the simulation without changing the computation power. The astrocyte agent is placed (or centered) within a grid space. The number of agents centered in any one grid space cannot exceed the maximum density parameter.

      At initialization, a random pairing of (x,y) is generated. If there is room (meaning that the maximum density parameter will not be exceeded, and that there are no fibers or microglial agents already in the grid space) to place an astrocyte agent in this grid space, one will be placed at (x,y). If there is no room, a new random pairing of (x,y) is generated and tested. This process continues until the number of astrocyte agents placed equals the initial astrocyte count.

States

      In the simulation, an astrocyte can be in one of four states:

      As soon as the concentration of IL-1B in the same grid space as the astrocyte agent exceeds the activation level, the agent's state changes from inactive to receptive. If the concentration ever falls underneath the activation level, the state of the astrocyte will return to inactive as long as it is not already blocking.

Absorption of IL-1B

      If the state of the astrocyte agent is not inactive, the agent can absorb IL-1B. Absorption is based on the receptor kinetics as described above. The Il-1B receptor binding rate is kf. It is used in conjunction with the IL-1B receptor equilibrium constant to determine kb. The IL-1B receptors per astrocyte is converted by the program (via a programmer defined conversion constant) into a concentration, R. Finally, the IL-1B receptor conversion rate defines k2. Using these parameters, the program can calculate k and h as defined in the receptor kinetics (case 2). Thus, if S is the concentration of IL-1B present, then the astrocyte will attempt to remove an amount, DS, on the macro time scale (time increment = DT) such that

DS = q*DT*k*S/(h+S)
where q is the astrocyte concentration. Having determined DS, the concentration of IL-1B in the same grid space as the astrocyte agent is decreased by DS while the concentration of IL-1B within the agent is increased by DS. Currently, there is no limit to how much IL-1B an astrocyte can absorb, although DS is limited to the amount of IL-1B present (an astrocyte cannot absorb more IL-1B than exists in the grid space).

Secretion of IL-6 and TNF

      Astrocytes secrete IL-6 and TNF based on the fraction of bound IL-1B receptors. As described in the section of neuron health, the fraction of bound IL-1B receptors is given by

Fraction = A/(h + A)

where A is the local concentration of IL-1B and h is derived from the full receptor model for the Michaelis-Menten kinetics. If the fraction of bound IL-1B receptors is greater than the fraction of bound IL-1B rec triggering IL-6 parameter, then the astrocyte will secrete IL-6 as described in the section on IL-6. Similarly, if the fraction of bound IL-1B receptors is greater than the fraction of bound IL-1B rec triggering TNF parameter, then the astrocyte will secrete TNF as decribed in the section on TNF.

Movement

      In order for astrocyte movement to occur, the state of the astrocyte must either be receptive or motile. In addition to having the proper state, the amount of IL-1B within the astrocyte agent must exceed the product of the IL-1B threshold for motility and the astrocyte concentration.

      The maximum speed parameter, v, is used to determine the number of macro time steps between movement, m. By letting DX be the length of a grid space and DT be the macro time increment, then m is the rounded integer value of (DX/v)/DT. Every astrocyte agent has an internal counter which keeps track of the number of macro time steps until movement can occur. Initially, the counter is set randomly between 0 and m (via a uniform distribution on the integer values). Every macro time step, the counter is decremented until it reaches zero. At this point, movement can occur and the counter is reset to m.

      Astrocytes try to move towards deposits of amyloid fiber. The astrocyte agent may move in one of several directions based on the grid. These directions are labeled 0-8 as follows

    0 1 2
    3 4 5
    6 7 8
where 4 determines the current position of the cell. Different rules apply to receptive and motile astrocyte agents.

      Receptive astrocytes sample the following grid spaces for amyloid fiber:

Motile astrocytes sample the following grid spaces for amyloid fiber:
Each numbered region corresponds to a direction of possible movement. The fiber concentrations in each region are normalized by dividing by the number of grid spaces sampled. The region of greatest fiber concentration is noted. The astrocyte will move in this direction in a Monte Carlo fashion with probability given by the probability of moving towards stimulus parameter. If the desired direction does not win, then a direction between 0 and 8 is chosen at random with a uniform distribution. If the astrocyte was receptive, its state will change to motile for the next time movement is to occur.

      Even if movement should take place, it can only occur provided that there is room for the astrocyte agent in the new grid space. Room in the new grid space is based on three criteria: (1) no microglia can be present, (2) no amyloid fibers can be present and (3) the number of astrocyte agents cannot exceed the maximum density parameter.

      If no fiber concentration is detected, then movement will not occur. If the state of the astrocyte is motile when this happens, then the state of the astrocyte will revert back to receptive.

Blocking

      Any astrocytes that are not inactive may become what we call blockers. In order for an astrocyte to become a blocker, some amyloid fiber must be in a neighboring grid space (these are the white spaces shown in the pictures for movement above). If fibers are present, then the astrocytes change the diffusivities in the immediate grid spaces normal to the direction at which the fiber was detected. The following pictures give examples of which grid spaces are affected given that the astrocyte agent is located at the center and the fiber is in the same grid space as the orange square.

         
For every whole astrocyte represented by an astrocyte agent, the chemical diffusivities of the affected grid spaces are reduced by a percentage equal to the reduction percentage for diffusivities parameter (call it r) with a probability equal to the probability of blocking parameter (call it p) in a Monte Carlo fashion. For example, let's say that the integer part of the astrocyte concentration parameter is 10. Then 10 uniformly distributed random variables between 0 and 1 are generated. If the value of the random variable is less than p, then we say that blocking for that astrocyte has "won" and the diffusivitiy needs to be reduced by r. If blocking "wins" 6 of the 10 times, then the new diffusivity, Dn, is calculated based on the old diffusivity, Do, such that
Dn = Do (1-r)6
for each of the three affected grid spaces.



Neurons

      Unlike microglia and astrocytes, neurons cannot move. The current formulation assumes that there is a neuron in every grid space. Earlier formulations assumed that multiple grid spaces lead to what was called a neuron block which is simply a population of neurons over some rectangular block of grid spaces. The computer code for this approach has not been deleted, but neither has it been updated.

      The major functions of neurons in the simulation are to absorb IL-1B, IL-6 and TNF, secrete soluble amyloid, and change health.

      At initialization, the neuron located in the grid space in the center of the environment begins with a concentration of IL-1B equal to the average of the maximum IL-1B absorbed parameter and the source triggering level parameter. In this manner, the center grid space will always be a source of soluble protein. This allows the feedback loop to begin and result in a meaningful simulation.

Absorption of IL-1B, IL-6 and TNF

      Absorption of IL-1B, IL-6 and TNF are all handled by the receptor kinetics described above. Therefore, for each chemical that is to be absorbed, there are corresponding user controlled parameters associated with them: the chemical receptor binding rate is kf; by multiplying it with the chemical receptor equilibrium constant, one obtains kb; the chemical receptor conversion rate is k2; and the chemical receptors per astrocyte is converted by the program (via a programmer defined conversion constant) into a concentration, R. Using these parameters, the program can calculate k and h as defined in the receptor kinetics (case 2). If S is the concentration of chemical (be it IL-1B, IL-6 or TNF) in the same grid space as the neuron, then the change in concentration, DS can be calculated on the macro time scale with time increment, DT, based on numerically solving the differential equation. Using a simple Euler method, we find that

DS = DT*k*S/(h + S).
DS is calculated for each chemical (IL-1B, IL-6 and TNF). The concentrations for each chemical are then updated such that the concentrations of the chemicals in the grid space are reduced and the concentrations within the neuron are increased by the appropriate DS.

      Special consideration is given when the chemical reaches the neuron's capacity to store it. The neuron will never absorb more chemical than it can hold. The limit for IL-1B is given by the maximum IL-1B absorbed parameter. The limit for IL-6 is given by the fatal IL-6 concentration parameter. The limit for TNF is given by the maximum TNF absorbed parameter. In the case where DS exceeds the room available to the chemical, DS is decreased to take the remaining room.

      In the unlikely case that DS exceeds S (the concentration of chemical that is present in the grid space), DS is changed so that it equals S. This is just a check to make sure that no negative concentrations of chemicals result from the absorption process. This makes sense both biologically and mathematically.

Secretion of Soluble Amyloid

      Neurons act as sources for soluble amyloid protein based on the amount of IL-1B they have absorbed. In order for a neuron to become a source of soluble amyloid, the concentration of IL-1B that it has absorbed must exceed the source triggering level parameter. At the time at which this occurs, a final test is done to see if it actually becomes a source based on the maximum proportion of neuron sources parameter, m. This "test" is done in a Monte Carlo fashion where a uniformly distributed random variable between 0 and 1 is generated. If its value is less than m, the "test" is passed and the neuron becomes a source of soluble amyloid. Because, IL-1B is not reduced within the neuron, this "test" is only done once and will determine whether the neuron will be a source for the duration of the simulation. Of course, neuron death causes any sources of soluble amyloid to die with the neuron.

      For details on the amounts of soluble amyloid secreted, see Amyloid Protein under Chemicals.

Health

      Neuron health varies as a rate depending on the proportion of bound receptors and local amyloid concentrations. Health varies between 0 and 1. A value of 0 for health indicates death. A value of 1 indicates that the neuron is in perfect health.

      Because neuron health varies as a rate, we describe changes in neuron health based on a differential equation. The first assumption we make is that as long as neurons do not decay past a certain point, the threshold of no recovery, they can recover via the logistic equation,

      We then assume that the cytokines, IL-1B, IL-6 and TNF can change health depending on the fraction of bound receptors. The fraction of bound receptors can be determined from the Michaelis-Menton kinetics. If A the concentration of chemical in the same grid space as the neuron, we can determine the fraction of bound receptors for the chemical, BA, as

where h is calculated the same way as was done for neuron absorption. In addition to using the fraction of bound receptors, we also define an effects on health parameter for each chemical, eA. We can combine these so that the new equation governing health is
where A denotes the chemicals, IL-1B, IL-6 and TNF. It is easy to see that if eA is negative, then the chemical is detrimental to health, while if it is positive, then the chemical is beneficial to health. If eA is zero, then that chemical has no effect on health.

      Finally we assume that the amyloid concentrations can affect health as well. We allow both the local concentration of soluble amyloid, S, and fibrous amyloid, F, to influence neuron health. Again we incorporate an effects on health parameter for soluble and fibrous amyloid, eS and eF, respectively. We make one final addition that scales the concentrations, so that the full effects won't be felt unless the concentration of S and F are at the amyloid conc affecting health, Smax, and fiber conc affecting health, Fmax, respectively. Combining all of these effects, the final equation for neuron health is

      We calculate neuron health on the macro time scale with time increment, DT, based on numerically solving the above differential equation, using a simple Euler method. H can never exceed 1 so that if the calculation is greater than 1, then H is automatically set to 1. Likewise, H can never be less than zero (the neuron is dead when H=0) so that if the calculation is less than zero, then H is automatically set to zero.

      Averaging the healths of all the neurons gives one the overall neuron health. Thus, if h(i) is the health of neuron i and there are n neurons, then the overall neuron health is simply





previous: Environment     next: Chemicals