Modeling spatio-temporal dynamics of network damage and network recovery

How networks endure damage is a central issue in neural network research. In this paper, we study the slow and fast dynamics of network damage and compare the results for two simple but very different models of recurrent and feed forward neural network. What we find is that a slower degree of network damage leads to a better chance of recovery in both types of network architecture. This is in accord with many experimental findings on the damage inflicted by strokes and by slowly growing tumors. Here, based on simulation results, we explain the seemingly paradoxical observation that disability caused by lesions, affecting large portions of tissue, may be less severe than the disability caused by smaller lesions, depending on the speed of lesion growth.


Introduction
The performance of networks in the case of damage is a central issue in biological as well as nonbiological networks [1,2].Intuitively the degree of disability following a damage must be correlated with the size of damage.However, different networks such as brain network, small world, scale free and random networks show different responses to spatial pattern of damage [3,4,5,6,7,8].Although there are many works ralated to network recovery [9,10,11], dependence of network recovery on network architecture in the case of spatial patterns of damage remains to be further explored.Another basic question is to consider the temporal patterns of damage.For example, a slow growing damage may lead to better recovery whereas a fast growing damage may lead to a recovery failure given a maximum size of spatial damage [12,13,14,15].There are experimental data on rat, cat and also monkey that shows a delay time between brain lesions has a strong effect on deficits that animals bear after operation [16,17,18,19,20,21,22].In a more clinical setting this proves to be a crucial aspect of how the brain recovers from strokes or some lesion growing tumors.The biological basis of recovery after stroke and how the brain reorganizes itself, for example, is still largely unknown [23,24] but certainly it is possible that the degree and speed of recovery varies considerably for different lesion locations and depends on structural alterations taking place in the spared brain tissue [16,17,18].Recently there has been some attempts using simple recurrent networks to show the effect of brain remapping and how it contributes to neuronal homeostasis following lesions [25,26].These findings, however, relate more to the spatial pattern of brain lesions and needs to be extended to what may be called a spatio-temporal pattern of brain injury.Here, a distinction has to be made between an acute stroke where there is a sudden damage inflicted to the brain and a slow growing damage such as a low-grade gliomas [12,13,14,15].Based on many experimental reports the slow growing injuries have a much better chance of being recovered than the injuries caused by acute lesions [16,17,18,19,20,21,22].As one might say if neurons can will, recovery comes as no surprise1 .In the present paper two simple networks, a three-layer and a homeostasis model with different temporal pattern of damages are studied.As will be shown a gradual injury is indeed resisted by these networks providing the lesions inflicted on the networks do not exceed a critical size.Our main finding is that a gradual injury will decrease the maximum amount of disability that the network may tolerate before it brakes down with no chance of recovery.

Three layer model 2.1 The model
In this model the nodes are binary with state 1 means active and state 0 inactive.The network learns the correct answer for each input.A measure of disability, the Hamming distance, is defined as the distance between the desired and the actual output.The network consists of three rings of nodes each containing N nodes.There is no connections between nodes in each layer but there is feed-forward connection from the previous to the next layer.Each node i is connected to a group of n neighbor nodes in the next layer (see figure1).

n N Input layer
Middle layer Output layer Figure 1: Three layer model network topology: each one of N nodes from input layer connects to n neighbor nodes in the process layer which are symmetrically asides the one which is in front of the node in the input layer.The same connections hold from process layer to the output layer.There isn't any connections between nodes of each layer but there is feed-forward connection from previous layer to the next one.
The fraction n N measures network completeness in the sense that if we remove a group of m neighbor nodes in the middle layer and m > n, the route from input layer to output layer for some nodes may be cut.For m < n, although the network may be initially disabled has the potential to recover.There is a sigmoid activation function for the middle layer and nodes in the output layer are binary threshold neurons.In each trial one input is shown to the first layer and the connection weights are changed according to a Hebbian-like learning rule as long as a difference exists between the desired and the actual output.
In the pre-lesion phase we start with random binary patterns in the input and output layer and give enough time for the system to learn.In the lesion phase we kill some nodes immediately or gradually in the middle layer and wait for recovery of the system.At each time step we calculate hamming distance in the output layer as a measure to quantify network disability.There are two approaches in the lesion phase.In the first one we remove the nodes gradually one by one until the total number of n l nodes are removed.In the other approach according to the resection experiments [16,17,18,19,20,21,22] the total number removed nodes, n l , is divided in to n p equal size packages where at each time step a package containing n l np nodes is removed.The simulation details and parameter values for this part can be find in the Appendix A.

Discussion on the learning rate, η
The learning rate was assume constant in 2. In a more realistic network η is held to be a function of network disability.We suppose that there is a threshold size of injuries above which the system may not be recovered.The η will decrease with the size of network disability as the injury may ruin the learning part of the system too.We suppose that the functionality of η with network disability behaves as in figure 2. In the homeostasis model described in section 3 the learning rate is also proportional to network disability.A more realistic relation between η and network disability.The η will decrease with the size of network disability as the injury may ruin the learning part of the system too.There is a threshold size that for injuries larger than that, the system could not be recovered anymore.

Three layer model with fixed η
The results for immediate and gradual injury for delay time "td" of 10 and total injury size of 80 is shown in figure 3a.In this diagram vertical axis indicates hamming distance as a parameter that shows the amount of network disability and horizontal axis is time.According to our three layer model, recovery/learning rate is a constant function of disability size.This is confirmed by the linear behavior of recovery shown in the red curve.In this diagram injury starts at time 1000.The red curve shows a sudden and large network disability which linearly recovers over time.In the blue curve although we start injury at time 1000,the nodes are removed gradually each 10 time step.Up to time 1200 the network somehow resists sudden huge disability but after 20 nodes being removed, the network suddenly shows a huge disability.From 1300 to 1800 when the last of 80 nodes is removed, there isn't a great change in network disability.At time 1800 the network starts to recover linearly and behave just like the red curve.There are two important things about this diagram.First, the maximum disability decreases as the nodes are gradually removed with increasing delay time.
There is however a new finding, to the effect that as we remove the nodes gradually the network disability doesn't increase gradually: Up to removing a specific number of nodes the network shows resistance but then suddenly changes its behavior and shows a great disability after which removing more nodes results in slower increase in disability.
This behavior points to a possible critical size of injury.To study that these results not just hold for a specific initialization we simulate the three layer model with the same parameters but for 200 different initializations.In the figure 3b the average of these simulations is shown.In this diagram the effects which we saw before is hold and also because of the small fluctuations of this diagram we observe new behavior in the recovery phase.

Three layer model with modified η
For better observation of this critical injury size effect we use the specific relation between learning rate η, and disability size described in section 2.2 (Also, see figure 2).With this dependency of η and hamming distance(disability size) we could capture the state of severe disability.In this study we set this critical hamming distance to 25.For hamming distances larger than 25 the learning rate η, is close to zero and therefore no learning is observed.With this newly defined η, we simulate three layer model for injury size of 30 with this newly defined η. Figure 3c show two single runs, one for immediate and one for delay time of 30 between removing nodes.In the first diagram the Maximum amount of disability (MAoD) never reaches 25 so the system could recover but in the next diagram where injury is immediate, At the begining the MAoD became larger than 25 so the recovery rate η drops down and system couldn't recover.In the figure 3d the same result for the average of 200 different initializations is shown.

Three layer model -Different kinds of injury
In the previous sections the behavior of three layer model for gradual injury has been shown briefly.
Here according to injury and recovery models four kind of sub-models has been developed: Type 1: removing each node gradually with a delay between removing them when we use a fixed η as a learning rate.Type 2: removing each node gradually with a delay between removing them when we use a modified η as a learning rate.Type 3: removing nodes in packages gradually with a delay between removing them when we use a fixed η as a learning rate.Type 4: removing nodes in packages gradually with a delay between removing them when we use a modified η as a learning rate.
In each sub-model the effect of delay has been studied.In the red curve where we remove the nodes suddenly the amount of network disability is more than in the blue one with gradual node removing.The existence of a critical injury size could be observed from the blue curve.During removing of the first 20 nodes the network shows a small disability (i).After which the hamming distance suddenly increases and the network shows a severe disability (ii).Removing more nodes from 1300 to 1800 doesn't change network disability a lot (iii).After removing last node, network start to recover linearly like the blue curve (iv).(b): All the the parameters in this diagram is equal to the figure 3a, but this one is the average of 200 different runs.The behavior in the lesion phase is the same as in figure 3a.(c): simulation with modified η for injury size of 30 for a single simulation.Red and blue curves are correspond to delay time (td) of 0 and 30 respectively.(d): simulation with modified η.For this curves we take the average of 200 runs with different initial conditions.The other parameters are the same as in figure 3c.

Type 1: No packages -Fixed η
For the maximum amount of 180 nodes and delays in the range of 0, sudden injury, to 1000 time steps, different systems had been simulated for 9 different initialization.A small subset of these simulations had been shown in red in the figure 4a.
For better demonstration of the results, by omitting the complex behavior of each system we briefly find the maximum amount of disability which plays the main role in our study and plot this feature in figure 5a.As it could seen, for large injury sizes and large delays the same MAoD "Maximum Amount of Injury" had been observed.In this section everything is like the previous section but the learning rate η which is we used a modified η as described in section 2.2.In the figure 4a the blue curves shows the behavior of this sub-model for some selected simulations.The same as previous section The MAoD is plotted in figure 5b.As it could seen, for large injury sizes and large delays the same MAoD "Maximum Amount of Disability" had been observed.

Type 3: packages -Fixed η
Based on the finding and experiments of previous works, [16,17,18,19,20,21,22] we know that removing the same amount of nodes in more than one package will decrease the MAoD too.To study this kind of injury in our model we set the injury size to 180 and remove them in different number of packages.For example if the number of packages is equal to four and delay set to 30, it means that every 30 time steps we remove one quarter of the whole injury size (45 nodes) together.
For the injury size of 180 nodes and delays in the range of 0, sudden injury, to 1000 time steps, different systems had been simulated for 9 different initialization.In the figure 4b a selection of this sub-model has been shown in red.The same as previous sections The MAoD is plotted in figure 5c.  3 Homeostasis model

The model
In search for a simple and more biological model to work on we select a model based on Butz [25,26].This model is based on placing simple spiking neurons on a ring where 80% of them are excitatory and the rest of them are inhibitory neurons.Neurons recive inputs on their dendrites and connect with other neuons through their axones.Here we have a total number of 100 neurons.For more information on network topology see [25,26].
Figure 6: Homeostasis model network topology: 80% of neurons (blue) are excitatory and the rest of them are inhibitory (red).Each neuron may has connections to its neighbors (Denderites and axones not shown).They can make synapse with nearer neighbors easier.
It has been shown [27,28] that for biological networks there is a moderate firing rate that may be different for each network.In this model we set this to 0.5 as we are not working on a specific network.The state of a network where the firing rate of neurons are similar to a moderate firing rate is called homeostasis.
There are two time steps in this model.One small time step called functional time step.In each functional time step we update the neuronal activity and the number of synapses dose not change.The second time step which is called morphological time step is the large time step.Here morphological time step is equal to 100 functional time steps.In each morphological time step according to neuronal activity in the last 1000 functional time steps we change the number of synapses in order to reach homeostasis.In each functional time step first we calculate X i the presynaptic input of each neuron i. X i simply is a summation over all excitatory and inhibitory neurons which connected to neuron i plus an external input which is from a poisson distribution.Then we update the firing probability F of each neuron with a sigmoid activation function.After updating F we update the state of each neuron.It may fire if a random number is smaller than F .In each morphological time step we change the number of synapses according to the history of F .The aim of the network is to reach homeostasis so we demand that the average of F over last 1000 functional time steps converges to 0.5.So if it was grater than 0.5 we change the number of synapses to decrease it.We do the opposite for dose neurons who have the average of F smaller than 0.5.For axones, first we calculate where ν is a small number that adjust the rate of converging to homeostasis, ∆F i = F i − 0.5 and A i is the number of connected axons of neuron i. changing in the number of dendrites is similar to previous relation.You can see the detailed relations for change in synaptic elements in [25,26].You can see in equation 1 that the changing rate is proportional to ∆F i .∆F i which shows how far is the neuron from the homeostasis can be interpreted as how bad this neuron works.By taking the variance of F over all neurons we have a parameter which shows how far is the whole network from homeostasis.The aim of this study is to see the relation between temporal patter of injury and network disability.In this model we suppose that the variance of F over all nodes in the past 1000 functional time step shows resembles the network disability.To study the effect of temporal pattern of injury, first we let the network to reach homeostasis then we start killing the neurons in the network by disconnecting them from the other network nodes.For the same number of nodes to remove, one time we remove them all at once and one time gradually and one by one with a delay time of "td" morphological time steps between each node removal.The simulation details and parameter values for this part can be find in the Appendix B.

Homeostasis model resuls
In figure 7 the result for an injury size of 10 is shown.In this figure the vertical axis shows the variance of F which demonstrate the network disability and the horizontal axis shows morphological time step.Injury starts at time 2700 where the network reaches homeostasis.For the blue curve all the 10 nodes are removed simultaneously but for red and green curves there are 10 and 20 morphological time steps respectively between removing each node.As clearly seen in the figure 7 the amount of network disability lowers as the delay time between removing nodes is increased.This result is in agreement with the results in the three layer model.
One important thing about homeostasis model is that the recovery rate is proportional to the size of disability.That is the reason all three different injury patterns recover almost at the same time.

Discussion
Using two simple network models we studied the effect of temporal pattern of node destruction on network recovery.For both networks increasing the time delay between removing the nodes lead to a decrease in the maximum amount of network disability.In the case of our simple three layer model the amount of network disability is much higher when nodes in the process layer are removed at once (figure 3).But when the node removal is more gradual the amount of disability is reduced.As shown in figure 7 the same results holds for the more biologically inspired model,the homeostasis model.In both cases big sudden injury results in serious disability but large delay between removing nodes do not show serious disability.In general these results could be accounted for as a generic property of a distributed network where in the case of a gradual removal of nodes the remaining nodes take part in the process of recovery.More complex interactions,however, may also be involved explaining the difference of a sudden versus gradual damage.We studied the effect of learning rate dependence on the size of damage and indeed it is reported that following damage there is a time period in which the process of recovery speeds up [14].Also, there are some animal models showing that a focal damage induces excitability and plasticity in the rest of the brain [29].Although the origins and mechanisms of these changes are still unknown it is worth studying these models in a more realistic setting, for example,it is known for a quite a while that brain recovery after sequential lesions depends on the amount of tissue resected at each surgical stage [30].In similar animal studies it is shown that in general a two-stage lesion has a much better chance of recovery than a one stage acute lesion [21].These are clear examples of what may be called a spatio-temporal pattern of brain damage.More importantly as has been argued by Daffauand et al. [14] the same sort of mechanism may be at work in aging up to a certain threshold where the system can no longer cope with the neural destruction.This indeed would be an interesting future research requiring a more realistic neural modeling.

A Simulation details for three layer model
In this simulation, One pattern which cover 50% of the input layer teach to the system and we demand a patter in the output which covers 50% of the output layer.For simulations with modified η, the value of η for hamming distance below h 0 is equal to η and for hamming distance above h 1 is rapidly decrease to η 0 (section 2.2).The activation function for the middle layer is where θ ml is threshold, β is noise and X t i = N j=1 w ij • j where w ij is the weight of connection from node j in the input layer to node i in the middle layer.Nodes in the output layer works on the basis of binary threshold model with threshold θ output .In this model the learining rule is Hebb-like and the change in connection wights ∆w is: where a, b and c are constants and η is the learning rate, s i and s j show the state of nodes i and j which is 1 if they fire and 0 otherwise.
In this study all the constants set as shown in table 1.

B Simulation details for homeostasis model
In this model the activation function is where θ is firing threshold and β is determine noise by changing the steepness of the sigmoid function.
In this study all the constants set as shown in table 2

Figure 2 :
Figure 2: A more realistic relation between η and network disability.The η will decrease with the size of network

Figure 3 :
Figure 3: (a):The amount of network disability (Hamming distance) is shown for sudden and gradual node destruction.In the red curve where we remove the nodes suddenly the amount of network disability is more than in the blue one with gradual node removing.The existence of a critical injury size could be observed from the blue curve.During removing of the first 20 nodes the network shows a small disability (i).After which the hamming distance suddenly increases and the network shows a severe disability (ii).Removing more nodes from 1300 to 1800 doesn't change network disability a lot (iii).After removing last node, network start to recover linearly like the blue curve (iv).(b): All the the parameters in this diagram is equal to the figure3a, but this one is the average of 200 different runs.The behavior in the lesion phase is the same as in figure3a.(c): simulation with modified η for injury size of 30 for a single simulation.Red and blue curves are correspond to delay time (td) of 0 and 30 respectively.(d): simulation with modified η.For this curves we take the average of 200 runs with different initial conditions.The other parameters are the same as in figure3c.

Figure 4 :Type 2 :
Figure4: (a): The effect of injury size and delay between removing the nodes had been shown.In each sub plot the amount of network disability as a function of time is shown for fixed (red) and modified η (blue).For larger delay times, even for large injury sizes the network shows very small disability.vertical grid line in each subplot shows time step 1000 which is the injury start time.(b): The effect of number of packages and delay between removing them had been shown.In each sub plot the amount of network disability as a function of time is shown for fixed (red) and modified η (blue).Increasing the number of packages and delay, decrease the MAoD.This diagram is an average of nine runs with different initial conditions.

Figure 5 :
Figure 5:The MAoD for different injury sizes and delays between removing nodes is plotted, using fixed (a) and modified η (b).The same is plotted for different number of packages and delays between removing them using fixed (c) and modified η (d).These diagrams are an average of nine runs with different initial conditions.Almost the same MAoD for different combinations of temporal and spatial pattern of injury.

Figure 7 :
Figure 7:The amount of network disability (shown by variance of F ) for different delay time between node removal is shown.By increasing the delay time "td" the maximum disability of network is decreased.Here blue line shows immediate defect but red and green ones show gradual defect with 10 and 20 time steps between node removal respectively.The maximum number of nodes removed in this diagram is 10.

Table 1 :
Data used in three layer model simulations

Table 2 :
[25,26]used in homeostasis model simulations You can see the more details for this model in[25,26].