Simulation of Spread and Control of Lesions in Brain

Abstract

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.

1. 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.

1.1. 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 defect-free immune system, in this world view, purportedly only attacks pathogens, the external agents that cause illness or disease [5–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–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.

2. Model

An undirected, fixed radius random graph 𝐺(𝑛,𝑟), with 𝑛 nodes (vertices) and radius of connectivity, 𝑟, 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 𝑟. 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 𝑑_𝑖 be the degree of the 𝑖th node, that is, the number of edges attached to it. The health status of each edge, at time 𝑡, is indicated by its “weight,” 𝑤(𝑗,𝑡) (𝑗=1,…,𝑑_𝑖), an integer number ranging from 0⋯𝑤_max. Edges with weight 𝑤_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 𝑗th edge is taken to be equal to 𝑤(𝑗,𝑡).

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 𝐩^𝑑_𝑖 (𝐩^𝑟_𝑖) of getting damaged (regenerated). In the general case, 𝐩^𝑑_𝑖 (𝐩^𝑟_𝑖) is a column vector of length 𝑤_max containing the transition probabilities from one state of health to another. Probability of programmed cell death, 𝑝_𝑝, is independent of the health status of the edge.

The functional or health status of the 𝑖th node is the sum over its edge weights, 𝑠_𝑖(𝑡)=∑^𝑑𝑖_𝑗=1𝑤(𝑗,𝑡). The maximum possible value of 𝑠_𝑖 is denoted by 𝑆_𝑖, which is realized when each 𝑤(𝑗,𝑡)=𝑤_max.

A node damaged by the pathological process generates an alarm signal when the ratio of its health status to the fully healthy state (𝑠_𝑖(𝑡)/𝑆_𝑖) falls below a threshold, 𝜏_𝑎𝑙. The signals received at the 𝑖th node are summed and propagated further when the summed signal strength reaches 𝑠_𝑖.

Programmed cell death is initiated at all the nodes where the propagated signals reach a threshold 𝜏_𝑏𝑓. The accumulated alarm signals in the region of programmed cell death, a circular region around the activated node of radius proportional to a parameter 𝐶_𝑏𝑓, 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 (𝑅_𝐼(𝑡)) 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, 𝛼×𝑅_𝐼(𝑡)×𝑅𝑂𝐼_𝑡=0, where 𝑅𝑂𝐼_𝑡=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 𝐶_𝑏𝑓×𝑅𝑂𝐼_𝑡=0.

3. Simulation

In the simulations reported here, a two-state model with 𝑤_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, 𝐩^𝑟_𝑖, 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 𝑛=400 and chosen a uniform random distribution of points in the unit square [0,1]×[0,1]. The radius of connectivity was set to 𝑟=0.2. All the results were also confirmed on a network of 𝑛=4000, with 𝑟=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 𝑡=0 in a region with 𝑅𝑂𝐼_𝑡=0=0.05 around the center at (0.5, 0.5); for 𝑛=4000, 𝑅𝑂𝐼_𝑡=0=0.015.

A uniform probability of pathologic damage 𝑝_𝑑=0.33 was used, with 𝛼=0.12. We varied 𝜏_𝑎𝑙, 𝜏_𝑏𝑓, and 𝐶_𝑏𝑓 to identify the conditions under which the pathological process could be controlled by the programmed cell death. Larger values of 𝜏_𝑏𝑓 indicate reduced sensitivity to the alarm signals whereas a larger value of 𝐶_𝑏𝑓 indicates that a larger area near the alerted node is subjected to programmed cell death. In the case of 𝜏_𝑎𝑙, larger values indicate quicker firing of alarm signals.

4. Results

Figure 1 shows the time series of damages caused to the system by both the pathological process and the programmed cell death process. The first column of panels in the figure shows the time course of instantaneous damages to the system. The middle column of panels shows the time course of the cumulative damages to the system. The last column of panels show the final state of the network at the end of the simulations.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 1 

Time course of damages to the system by the pathological and programmed cell death processes. The first column of panels shows, separately, the instantaneous damages due to both the processes; dotted lines with asterisks indicate the damages due to the pathological process, and bold lines with filled circles indicate damages due to programmed cell death. The second column of panels shows the cumulative damages with time. Again, the damages effected through both the processes have been separately shown (same symbols as earlier), as also the sum total damages to the system (square symbol). The last column of panels (color online) show the state of the network at the end of the simulations; the dotted lines indicate healthy edges (axons), full (blue) lines indicate edges damaged due to programmed cell death, and dark (red) lines indicate edges damaged due to the pathological process. For all the panels, 𝜏𝑏𝑓=0.5 and 𝜏𝑎𝑙=0.7, while the 𝐶𝑏𝑓 values, for each row, top to bottom, are 0.2, 0.8, and 1.5, respectively.

There are three typical scenarios which are illustrated in Figure 1, in the three rows from top to bottom. Figures 1(a)–1(c) show a scenario where the programmed cell death is not of sufficient strength to significantly affect the pathological process. Note that the instantaneous damages from 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 𝜏_a𝑙, high 𝜏_𝑏𝑓, and low 𝐶_𝑏𝑓 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 𝐩^𝑑_𝑖 for that edge. Secondly, the random network itself is generated by the random distribution of the 𝑛 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.

(a)

(b)

(c)

(d)

Figure 2 

Phase portraits of the dynamics. The 𝑦-axes in Figures 2(a)-2(b) and the gray (color online) scale in Figures 2(c)-2(d) indicate the fraction of the total number of simulations where the pathological process was successfully arrested. In Figure 2a, 𝜏𝑎𝑙 has been held fixed at the shown value while the different curves are, from left to right, for 𝜏𝑏𝑓=0.1,0.3,0.5,0.7, and 0.9. In Figure 2(b), the different curves, from right to left, are for 𝜏𝑎𝑙=0.9,0.7,0.5, 0.3, and 0.1; 𝜏𝑏𝑓=0.5 is held fixed. The 𝑆-curves that result with increasing 𝐶𝑏𝑓 values indicate the smooth nature of the transition in dynamical behavior in the parameter space.

From Figure 2, we see that the transition from uncontrolled pathological process to arrested pathological process is smooth as 𝐶_𝑏𝑓 is varied from low to high values. In Figure 2(a), 𝜏_𝑎𝑙 has been held fixed and the different curves, from left to right, are for different 𝜏_𝑏𝑓 values, from 0.1 to 0.9 in steps of 0.2. In Figure 2(b), 𝜏_𝑏𝑓 has been held fixed and the different curves are, from right to left, for 𝜏_𝑎𝑙 values ranging from 0.1 to 0.9, in steps of 0.2. We shall denote by 𝐶^𝑐𝑟_𝑏𝑓(𝜏_𝑎𝑙,𝜏_𝑏𝑓) the critical value of 𝐶_𝑏𝑓 at which these 𝑆-curves attain a value of 1, that is, all instances of simulations result in the growth of pathological process being arrested. The combined picture in the three parameter space is presented in Figures 2(c)-2(d), from two different perspectives. The three different scenarios presented, from top to bottom, in Figure 1 indicate, respectively, the three different parts, from left to right, of a typical 𝑆-curve of Figures 2(a)-2(b).

As seen from Figure 2, pathological process is always controlled if 𝐶_𝑏𝑓>𝐶^𝑐𝑟_𝑏𝑓. Nevertheless, the sum total damage to the system is not the same for all values of 𝐶_𝑏𝑓>𝐶^𝑐𝑟_𝑏𝑓; in fact, the damage is greater, the larger the value of 𝐶_𝑏𝑓. 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 𝐶_𝑏𝑓 values have been plotted in Figure 3. Three different curves for three different 𝜏_𝑏𝑓 values are shown in this figure; similar graphs can be constructed for different 𝜏_𝑎𝑙 values as well (not shown here; see [16]). These averages have been taken at 𝑡=20 in each case. For 𝐶_𝑏𝑓<𝐶^𝑐𝑟_𝑏𝑓 values, the damages due to the pathological process as well as sum total damage are still growing and have not become stationary at 𝑡=20; for 𝐶_𝑏𝑓>𝐶^𝑐𝑟_𝑏𝑓 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 𝐶_𝑏𝑓=𝐶^𝑐𝑟_𝑏𝑓. For 𝐶_𝑏𝑓>𝐶^𝑐𝑟_𝑏𝑓, the programmed cell death is clearly effecting more damage than is necessary to arrest the pathological process. Since 𝐶^𝑐𝑟_𝑏𝑓 depends on 𝜏_𝑏𝑓 and 𝜏_𝑎𝑙, it is not surprising to see that (cf. Figures 2(a) and (3)) lesser damage results when 𝜏_𝑏𝑓 is small.

Figure 3 

Optimality in sum total damage to the system while effecting arrest of pathological process is shown here. With increasing 𝐶𝑏𝑓 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 𝐶𝑏𝑓 with larger 𝜏𝑏𝑓, and the minimum value also shifts upward. A similar situation occurs with different 𝜏𝑎𝑙 values (not shown here). Note that the 𝑦-axis values have been taken at 𝑡=20.

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, 𝛼×𝑅_𝐼(𝑡)×𝑅𝑂𝐼_𝑡=0. This is achieved in all instances of simulation when 𝐶_𝑏𝑓>𝐶^𝑐𝑟_𝑏𝑓. Currently, mathematical analysis of this feature is being carried out to establish the relationship of 𝐶^𝑐𝑟_𝑏𝑓 with 𝜏_𝑎𝑙 and 𝜏_𝑏𝑓, and the results will be reported soon.

5. 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.

Acknowledgments

This work was carried out during a period spent at the Department of Pharmaceutical Sciences, State University of New York at Buffalo (UB), USA. Encouragement, discussions, and funding from Dr. Murali Ramanathan are acknowledged. Use of the computational facilities at the Center for Computational Research, UB, facilitated the work.