Simulation of Spread and Control of Lesions in Brain

A simulation model for the spread and control of lesions in the brain is constructed using a planar network (graph) representation for the central nervous system (CNS). The model is inspired by the lesion structures observed in the case of multiple sclerosis (MS), a chronic disease of the CNS. The initial lesion site is at the center of a unit square and spreads outwards based on the success rate in damaging edges (axons) of the network. The damaged edges send out alarm signals which, at appropriate intensity levels, generate programmed cell death. Depending on the extent and timing of the programmed cell death, the lesion may get controlled or aggravated akin to the control of wild fires by burning of peripheral vegetation. The parameter phase space of the model shows smooth transition from uncontrolled situation to controlled situation. The simulations show that the model is capable of generating a wide variety of lesion growth and arrest scenarios.


Introduction
MS affects about one million people worldwide and causes physical and cognitive disability. There are three types of MS, relapsing-remitting, secondary progressive, and primary progressive, that differ in the dynamical patterns of disease progression. There are as yet no known cures for MS. Patients with relapsing MS are currently treated with drugs that exert immunomodulatory effects and slow the progression of the disease; there are no effective treatment options for the progressive forms of MS [1,2].
MS is postulated to be a cell-mediated autoimmune disease directed against myelin components of the CNS. Myelin is an electrically insulating phospholipid layer that surrounds the axons of many neurons. The disease is characterized by both inflammatory immune responses and neurodegeneration. The prevailing hypothesis on MS pathogenesis is that autoreactive T-lymphocytes, a cell type in the immune system, orchestrate a complex cascade of events that cause blood-brain barrier disruption and invasion of immunologically aggressive cells into the CNS. However, the exact causes of MS still remain unknown [3,4]. The long-term goals of this research are to develop disease models that can be used to evaluate therapeutic strategies for this disease and, in this report, the specific focus is on evaluating a network model for MS lesion dynamics. Literature survey indicates that network approaches have not been studied extensively for disease modeling in MS.

Previous Work.
Conventional models for autoimmunity are premised on the occurrence of defects in the immune system that cause it to turn against the host tissue. A defectfree immune system, in this world view, purportedly only attacks pathogens, the external agents that cause illness or disease [5][6][7]. However, an alternative viewpoint has been advocated where auto-immunity is seen as the usual immune response, but directed against those components of the body which, in normal conditions, are inaccessible to the immune system [8][9][10][11][12][13][14]. For example, in the danger model, developed by Matzinger [10,11], it is posited that stressed and injured tissues can mediate immune responses through the generation of appropriate "danger" signals. This is as opposed to the activation through recognition of external pathogenic cell types from host tissue in the conventional models. The concept of comprehensive immunity, developed by Nevo et al. [12,13] complements this alternate perspective; experimental results supporting their idea have also been reported [14]. The present network model is inspired by the alternative viewpoint. The key elements of the model consist of a pathological process that causes cellular damage and programmed cell death (apoptosis) initiated through an intercellular signaling component. The programmed cell death deprives the pathological process of healthy tissue which is necessary for its propagation in space and time. In this, it resembles the action of firemen who burn peripheral vegetation to contain forest fires. Inter-cellular signaling is a key feature of the model that allow pathologically damaged cells to propagate alarm signals and initiate programmed cell death.

Model
An undirected, fixed radius random graph G(n, r), with n nodes (vertices) and radius of connectivity, r, is constructed to represent the CNS in this 2D network model. Fixed radius implies that nodes are connected only if they are within a distance of r. Biologically, the nodes of the graph can be viewed as representing cell bodies or functional units and the edges (bonds) of the graph can be viewed as axons or the interconnections between functional units.
Let d i be the degree of the ith node, that is, the number of edges attached to it. The health status of each edge, at time t, is indicated by its "weight," w( j, t)( j = 1, . . . , d i ), an integer number ranging from 0 · · · w max . Edges with weight w max are fully functional or healthy units (as at the beginning of simulation), and those with weight zero, are dead. Extending the same logic, the amplitude of the signal propagated along jth edge is taken to be equal to w( j, t).
In the pathological process, the edges are damaged by lowering their weight by a single unit. However, in the programmed cell death process, edge weights are directly reduced to zero. In the regeneration process, edge weights are raised by a unit.
The pathological and regeneration processes are driven by probabilistic events wherein each edge in the affected region, in each time unit, has a certain probability p d i (p r i ) of getting damaged (regenerated). In the general case, p d i (p r i ) is a column vector of length w max containing the transition probabilities from one state of health to another. Probability of programmed cell death, p p , is independent of the health status of the edge.
The functional or health status of the ith node is the sum over its edge weights, A node damaged by the pathological process generates an alarm signal when the ratio of its health status to the fully healthy state (s i (t)/S i ) falls below a threshold, τ al . The signals received at the ith node are summed and propagated further when the summed signal strength reaches s i .
Programmed cell death is initiated at all the nodes where the propagated signals reach a threshold τ b f . The accumulated alarm signals in the region of programmed cell death, a circular region around the activated node of radius proportional to a parameter C b f , get reset to zero. No additional signals are generated at these nodes to the alarm signals generated in the pathologic process.
The spread of the pathologic process is driven by the success rate in causing cellular damage. The fraction of edges (R I (t)) damaged in a particular time step, among the total number of healthy edges visited, is the rate of damage due to the pathologic process. The rates of damage due to the pathologic and the programmed cell death are computed in terms of the initial lesion size so that the final results are invariant with respect to the initial lesion size. Thus, the radius of the region affected by the pathologic process increased or decreased according to the formula, α×R I (t)×ROI t=0 , where ROI t=0 is the radius of the region at the center where the initial lesion is seeded. In a similar fashion, the region of programmed cell death was computed as C b f × ROI t=0 .

Simulation
In the simulations reported here, a two-state model with w max = 1 has been employed, that is, there are no intermediate states of health, and the edges are either alive or dead. Additionally, the regeneration probability, p r i , was set to zero in order to focus exclusively on the effects of the interplay between the pathological and programmed cell death processes on lesion structure and dynamics. A few preliminary results using such a configuration was reported earlier [15].
We have set n = 400 and chosen a uniform random distribution of points in the unit square [0, 1] × [0, 1]. The radius of connectivity was set to r = 0.2. All the results were also confirmed on a network of n = 4000, with r = 0.06. Average degree strengths of the order of 10 are obtained in these configurations; degree distribution is Gaussian. The pathological process was initiated at t = 0 in a region with ROI t=0 = 0.05 around the center at (0.5, 0.5); for n = 4000, ROI t=0 = 0.015.
A uniform probability of pathologic damage p d = 0.33 was used, with α = 0.12. We varied τ al , τ b f , and C b f to identify the conditions under which the pathological process could be controlled by the programmed cell death. Larger values of τ b f indicate reduced sensitivity to the alarm signals whereas a larger value of C b f indicates that a larger area near the alerted node is subjected to programmed cell death. In the case of τ al , larger values indicate quicker firing of alarm signals.  programmed cell death are hardly ever above zero. Also, it is seen from Figure 1(b) that the contribution of programmed cell death to the sum total of damages is insignificant. This situation occurs with a suitable combination of low τ al , high τ b f , and low C b f values. Figures 1(d)-1(f) show a slightly more complex situation. In this case, programmed cell death is clearly the dominant effect. The instantaneous damages caused by both the processes are consistently nonzero (Figure 1(d)) and the cumulative damages (Figure 1(e)) continue to grow. The total damage, thus, continues to spread. In Figures 1(g)-1(i), the pathological process has been well controlled. The instantaneous damages have fallen to zero in Figure 1(g), and the cumulative damages (Figure 1(h)) have leveled off. The final state of the network (Figure 1(i)) shows that the damage is also minimal in terms of the fraction of edges damaged. As seen from Figures 1(a), 1(d), and 1(g), the time series is stochastic. There are essentially two sources of randomness in the model. Firstly, the pathological process is simulated by a binomial process wherein each edge visit will lead to successful damage if the generated random number falls below the value in p d i for that edge. Secondly, the random network itself is generated by the random distribution of the n points in the plane. The complete picture of the transition from uncontrolled growth of the pathological process to the situation where the pathological process has been well arrested is seen in the parameter phase space graphs shown in Figure 2, where an averaging has been effected over the two sources of randomness. The phase space diagrams are the results of averaging over ten different networks, with the dynamics averaged over a thousand iterations.

Results
From Figure 2, we see that the transition from uncontrolled pathological process to arrested pathological process is smooth as C b f is varied from low to high values. In Figure 2(a), τ al has been held fixed and the different curves, from left to right, are for different τ b f values, from 0.1 to 0.9 in steps of 0.2. In Figure 2(b), τ b f has been held fixed and the different curves are, from right to left, for τ al values ranging from 0.1 to 0.9, in steps of 0.2. We shall denote by As seen from Figure 2, pathological process is always controlled if C b f > C cr b f . Nevertheless, the sum total damage to the system is not the same for all values of C b f > C cr b f ; in fact, the damage is greater, the larger the value of C b f . Clearly, it is desirable to effect control of the pathological process with the least sum total damage to the system. With this in mind, average fractional damages at different C b f values have been plotted in Figure 3. Three different curves for three different τ b f values are shown in this figure; similar graphs can be constructed for different τ al values as well (not shown here; see [16]). These averages have been taken at t = 20 in each case. For C b f < C cr b f values, the damages due to the pathological process as well as sum total damage Average fractional damage Figure 3: Optimality in sum total damage to the system while effecting arrest of pathological process is shown here. With increasing C b f values the minimum in sum total damage to the system occurs at about the same value as the critical value at which the fraction of simulations in which the pathological process is arrested attains unity (cf. Figure 2(a)). This minimum occurs at higher values of C b f with larger τ b f , and the minimum value also shifts upward. A similar situation occurs with different τ al values (not shown here). Note that the y-axis values have been taken at t = 20.
are still growing and have not become stationary at t = 20; for C b f > C cr b f values, the averages have become stationary. Nevertheless, these curves indicate that the least sum total damage to the system, with pathological process arrested, is obtained at C b f = C cr b f . For C b f > C cr b f , the programmed cell death is clearly effecting more damage than is necessary to arrest the pathological process. Since C cr b f depends on τ b f and τ al , it is not surprising to see that (cf. Figures 2(a) and (3)) lesser damage results when τ b f is small.
From the above, it is clear (cf. Figures 1(d)-1(e)) that arrest of the pathological process does not necessarily occur if the damage due to the programmed cell death process is greater than the pathological process. What is necessary [16] is that the programmed cell death process be able to encircle the region affected by the pathological process, and, furthermore, be able to create an envelope region of sufficient thickness to offset its likely growth factor, α × R I (t) × ROI t=0 . This is achieved in all instances of simulation when C b f > C cr b f . Currently, mathematical analysis of this feature is being carried out to establish the relationship of C cr b f with τ al and τ b f , and the results will be reported soon.

Conclusions
A physically motivated 2D network model was developed for the CNS and employed to study the process of lesion formation and spread in MS. Intercellular signalling of distress by the damaged cells is a key feature of the model which leads to programmed cell death getting activated in an attempt to arrest the lesion progress. The model demonstrates that the spread of the pathologic process can be arrested by programmed cell death when the geometry of the damage inflicted by the latter leads to an envelope, of sufficient thickness, being created encircling the area of pathological process. Such an envelope of dead cells deprives the pathological process of healthy cells which can sustain its growth. The model shows a smooth transition, as parameters are varied, from the situations of run-away pathological process, through aggravated damage to the system caused by unsuccessful firing of programmed cell death, to the creation of successful envelope around the pathological process.
The model complements the alternate viewpoint on autoimmunity which posits that cells and tissues signal distress and activate the immune system. Such a viewpoint circumvents the need for the immune system to store information about likely pathogens and, also, makes it capable of acting in instances of cellular damage resulting from nonpathogenic causes. Further study of the model along with identification of the possible biological constituents should enable comparisons with experiments and a more detailed exposition.