Function Based Brain Modeling and Simulation of an Ischemic Region in Post-Stroke Patients using the Bidomain

BACKGROUND
Several studies have shown that post-stroke patients develop divergent activity in the sensorimotor areas of the affected hemisphere of the brain compared to healthy people during motor tasks. Proper mathematical models will help us understand this activity and clarify the associated underlying mechanisms. New Method. This research describes an anatomically based brain computer model in post-stroke patients. We simulate an ischemic region for arm motion using the bidomain approach. Two scenarios are considered: a healthy subject and a post-stroke patient with motion impairment. Next, we limit the volume of propagation considering only the sensorimotor area of the brain. Comparison with existing methods. In comparison to existing methods, we combine the use of the bidomain for modeling the propagation of the electrical activity across the brain volume with functional information to limit the volume of propagation and the position of the expected stimuli, given a specific task. Whereas just using the bidomain without limiting the functional volume, propagates the electrical activity into non-expected areas.


RESULTS
To validate the simulation, we compare the activity with patient measurements using functional near-infrared spectroscopy during arm motion (n=5) against controls (n=3). The results are consistent with empirical measurements and previous research and show that there is a disparity between position and number of spikes in post-stroke patients in contrast to healthy subjects.


CONCLUSIONS
These results hold promise in improving the understanding of brain deterioration in stroke patients and the re-arrangement of brain networks. Furthermore, shows the use of functionality based brain modeling.


Introduction
According to the World Health Organization (WHO), stroke is the second leading cause of death and the third leading cause of disability worldwide (Organization et al., 2000(Organization et al., -2012. Stroke may lead to motion impairment in adults, such as post-stroke hemiplegia or hemiparesis caused by the effects of focal lesions in the brain Rathore et al. (2001). When a stroke occurs, a portion of the brain could be affected and impair motor functions, generally within one hemisphere Bonita and Beaglehole (1988). Most stroke survivors present motion impairment in their arms and there is high variability in the recovery process. Hence, neuromuscular modeling and simulation have vast potential to improve patient care by helping to elucidate the cause-effect relationships in subjects with neurological and musculo-skeletal impairments and establish effective rehabilitation treatments Cramer (2008).
The early mathematical models of the brain, modeled the spatiotemporal dynamics of the brain's electrical activity, simulating sources as dipoles following Poisson's equation He et al. (2002), He (2010). More recent models, focus on the relationship between the neuron interactions governed by ion currents and stimuli at a microscopic scale and the brain as a whole at a macroscopic scale considering reactiondiffusion systems Somogyvári and Érdi (2019). One of these models, is the one described in the Virtual Brain.
The Virtual Brain is a complex network characterized by functional structural connections, and it is able to simulate the brain's activity at different scales Sanz Leon et al. (2013). For example, one previous study Falcon et al. (2015) used the Virtual Brain to simulate brain activity in stroke patients to determine potential therapeutic interventions. The model in the Virtual Brain depends on structural connections, but after brain damage, there may be altered conductivity of the ischemic region Abboud et al. (1995), changing the connections. Thus, as an alternative, we propose using the bidomain model to simulate electrical activity in the brain.
The bidomain model assumes that the electrical activity in excitable tissue is generated by the depolarization of the cell membrane between the intracellular and extracellular domains, where both domains share the same space Sundnes et al. (2007), Henriquez (1993). The bidomain approach was developed to describe electrical activity in cardiac tissue; but since then, it has been adapted to other systems like the brain Szmurlo et al. (2006Szmurlo et al. ( , 2007, Yin et al. (2013),  and skeletal muscles Pascual-Marqui (1999), Weinstein et al. (1999) and . This model has the advantage of combining the overall electrical activity in an organ and microscopic nonlinear cell model ion interactions like the current-voltage dynamics given by cell models, such as the Fitzhugh-Nagumo Izhikevich and FitzHugh (2006) and Hodgkin-Huxley Noble (1962) models. In contrast, with the connections representation in other models, the activity diffusion is given by fiber directions. Each point in the mesh solves a node-wise cell model; hence, through a diffusion tensor calculates the overall electrical activity in a two-variable system of partial differential equations (PDEs). In addition, the simplified version of the bidomain model, the monodomain model Potse et al. (2006), reduces the system from a two-variable PDE system to one Sundnes et al. (2006). This simplification, allows us to use the finite element method (FEM) to require less CPU-power to simulate a volume section of the brain Dumont et al. (2019), and the effect of a ischemic region.
Human motion, whether voluntary or involuntary, is developed by spatial and temporal patterns of muscle contractions, controlled by neural circuits in the brain and spinal cord. Upper motor neurons in the primary motor cortex control voluntary motion Rizzolatti and Luppino (2001). E.g., in arm motion, the corresponding area of the contra-lateral hemisphere generates electrical signals to execute movement. Then, these signals propagate into the brain stem and spinal cord. Next, these signals descend through the corticospinal tract in the motor pathways until reaching the corresponding lower motor unit Purves et al. (2008). In post-stroke patients, the upper motor neurons in the primary motor cortex activation is changed, often preventing patients from moving their limbs Shao et al. (2009).
This study develops the process of simulating motion impairment in post-stroke patients by modeling a brain injury. We simulated a brain ischemic region by modifying the electrical conductivity of some mesh elements, similar to what is reported in Abboud et al. Abboud et al. (1995), but keeping the fiber directions organization.
For validation, we compare the global simulated activity in the brain between 5 post-stroke patients and 3 healthy subjects executing a motor task and measuring the brain activity using functional near-infrared spectroscopy (fNIRS). fNIRS is a specialized research imaging technique that uses near-infrared light to measure changes in the hemodynamic characteristics of the brain and allows to register concentration changes in deoxy-hemoglobin, oxy-hemoglobin, and total       A. Lopez -Rincon, et al. Journal of Neuroscience Methods 331 (2020) 108464 hemoglobin (HbT) Huppert et al. (2006). For our analysis, we used the HbT signal, as it is considered to be proportional to cerebral blood volume and brain activation Krieger et al. (2012) and Jasdzewski et al. (2003).

Bidomain Formulation for Electrical Activity Propagation
The bidomain model is described by the following equations Tung (1978): (1) where v is the transmembrane potential, u e is the extracellular potential, M i is the intracellular conductivity tensor, and M e is the extracellular conductivity tensor. M i is a 3x3 size diagonal matrix with diagonal values σ l , σ t , and σ n . C m is the transmembrane capacitance, β is the membrane surface-to-volume ratio, I ion is the total ionic current, I app is the external applied current source, and w represents the ionic model variables. For our simulations, we use the Fenton-Karma model with three variables Fenton and Karma (1998). If we assume a linear relationship between the intra and extra-cellular conductivity tensors such as M e = λM i , where λ is a constant scalar, we can reduce the bidomain model to the monodomain: with the boundary condition: and the following equation to get the extra-cellular potential: We then scale the equations with: in order to simplify the notation.

Modeling the Brain and Stimuli
The bidomain and monodomain models depend on the fiber directions for approximating the propagation of electrical activity. To create the fiber directions in the brain, we consider the white matter and draw its fiber directions tangent to the surface by taking into account its anisotropy (the gray matter anisotropy is negligible, Hallez et al. (2007)); see the fiber directions in Fig. 1.
The execution of voluntary movements begins in the motor cortex, which is the outermost layer. Hence, we focus our stimuli there, which is analogous to I app in Eq. (3). The position of the stimuli comes from the selection of the points given by hand motion and execution of motion in the LinkRBrain database (Fig. 2). LinkRBrain contains the information of thousands of peaks and coordinates from blood-oxygen-level dependent signals, obtained by functional magnetic resonance imaging, extracted from several published neuroimaging studies Mesmoudi et al. (2015). The stimuli is shown as the central point of the maroon volume shown in Fig. 3.

Simulating Functional Section
The bidomain and monodomain models are useful in tissue cells analysis via gap junctions, considering them as a syncytium. A classical   A. Lopez -Rincon, et al. Journal of Neuroscience Methods 331 (2020) 108464 example of this is the heart Sundnes et al. (2007). In contrast, in the brain neural signals propagate from one part of the brain to another, through discrete axons while neighboring axons and neurons remain at rest. In the case of motor movement, only a fraction of the brain is activated in comparison to the whole volume. Thus, to improve the simulation, we reduce the simulated volume to the sensorimotor cortex Jack et al. (1994), and to the left hand movement using LinkRbran Mesmoudi et al. (2015) (Fig. 4).

Ischemic Region
Ischemic tissue has very low conductivity Clay and Ferree (2002). Therefore, to model a focal stroke lesion, we modify the properties in certain elements from the brain mesh. In practice, we reduce the conductivity tensor values at several elements in the mesh as shown in Fig. 3, where the green volume represents the ischemic region. A similar approach is used in Swanson et al. to simulate a glioma Swanson et al. (2002).

Numerical Example
To test our setup, we reduce the brain system considering I ion and I app as the source function f(x, t) with σ l = σ t = σ n = 1 and = + a . Thus, reducing the system to: where g(x, t) and h(x, t) are the Dirichlet and Neumann boundary conditions; this is equivalent to the heat equation with a source. To test our system we use two spheres with a mesh of 17056 elements (Fig. 5). We divide the boundary, being the inner sphere ∂B 1 Dirichlet conditions and the outer sphere ∂B 2 Neumann conditions. The analytical solution is: for t 0 < t < t max . The error is given by the norm v || || 2 between the analytical solution and the calculated values at t max . We run the system from t 0 = 1 to t max = 2.5 with Δt = 0.01. The comparison of the calculated and the analytical isovalues solution at t max is shown in Fig. 6 for an error at t max of 0.00536. Fig. 15. Comparison between the measurements in the control subject (top) and the post-stroke patients for the different tasks. Column 1; assisted movement column 2; movement of both arms and column 3; self-assisted movement of the arm. A. Lopez-Rincon, et al. Journal of Neuroscience Methods 331 (2020)

Non-Functional Brain Simulation Comparison
Given the same parameters, stimuli, and fiber directions (Fig. 1), we simulate the electrical activity in the brain for the healthy subject and post-stroke patient with a focal lesion. In this case we consider the whole brain volume as a syncytium, and as expected the electrical activity propagates to functional areas, beyond what we expected; see Fig. 7.

Functional Brain Simulation
To make a more realistic approach to the brain's electrical propagation, we reduce the volume propagation to a functional volume, based in biological behaviour. Thus, we simulate the electrical activity, but we limit the volume of propagation as explained in section 2.3, to only the sensorimotor cortex; see Figs. 8 and 9 .

Comparison with Measurements
In order to verify the above simulations, we measured 5 post-stroke patients brain activity and 3 healthy controls with fNIRS. The representative position of the fNIRS leads is shown in Fig. 10.
The patients were asked to perform 3 tasks movements: First, to raise both arms (BA) Fig. 11. Next, to raise both arms with therapist assistance for the affected arm (AA) Fig. 12. Finally, to raise both arms using self-assisted arm movement, where the subject raises the affected arm with the non-affected arm (SA) Fig. 13. The control subject performed the movement with the dominant arm. The protocol consisted of 20 sec. rest plus 20 sec. for the task with 10 repetitions each.
Next, we compared the sum of the activity of the HbT signal for the different tasks in the patients Fig. 14, and controls Fig. 15. Finally, we compared the measured activity with the overall activity in our simulations (Figs. 16 and 17 ). In addition, we compared the projected electrical activity into the whole brain volume in Fig. 18.
To compare the simulated and measured activity in the Fig. 14-17, we remove the signals offset and consider the number of peaks in each hemisphere for patients (Table 1 ) and controls (Table 2). For our simulation, we simulate only in one hemisphere Fig. 18. The simulated healthy brain presents 1 peak, and the ischemic simulation 2 peaks.

Discussion
The correlation between the measured hemodynamic response and simulated membrane voltage is given by a spatial map that correlates  the excitatory membrane activity and the HbT signal, such as the one presented by Ma et al. (2009), which relates the balloon model of the hemodynamic response and the electrophysiological model. As expected, there is a divergence in post-stroke patient simulation and measurements with controls.
The simulated results are consistent with measurements from the literature and extensive studies. For example, Mihara and Miyai (2016) noted that patients with a post-stroke condition and chronic degenerative ataxia showed higher prefrontal activation in a steady gait study. A similar association was observed by Takeda et al. (2014), who found patients with moderate hemiparesis present increased activation during affected-hand grasping or higher intra-hemispheric activity in the brain Baldassarre et al. (2016).
From the visual inspection of Figs. 14-17, it is clear that the activity is concentrated for the control subject, whereas the activity expands for post-stroke patients. This is consistent with the findings of Falcon et al. (2015), who identified re-arrangement in the brain networks after a stroke. As a comparison, the HbT signal is consistent with the measurements from the dataset of Buccino et al. (2016), where a healthy subject is asked to move right and left arms and is measured with a 34 channels fNIRS; see  show there is one source for the simulated healthy subject, whereas there are two for the simulated post-stroke ones. This relationship is similar to the results between the control and the patients in Figs. 14, 15. A numerical explanation obtained from simulations is that the energy builds up around the focal lesion, which in turn makes close elements to depolarize in different directions. Therefore, generating more peaks in total, even if the only change between the two systems is the conductivity decrease in some mesh elements.
This relationship is more clear when comparing the number of peaks between control 1 and patient 2 Fig. 20, with the simulation Fig. 21. In both cases an extra peak appears.

Conclusion
In this research we propose a computational model to simulate the effects of focal brain lesions in post-stroke patients for hand motion. As mentioned before, this is an abstraction of the inner-workings of the brain, given we do not consider the subnetworks in the brain that are involved in motion and as a result, it is a reduced system that allows us to perform the simulation.
This method uses a numerical simulation with FEM in an open loop control model, which simulates the brain activity with and without a focal lesion. The simulations and comparisons with real measured data show that this model could be a suitable approach to modeling human motion and blockages in the brain originated from strokes or chronic diseases. In addition, it is a step forward to build a complex electrophysiological simulator in order to state functions of brain activity and thus, develop new therapies.

Implementation
Several applications were created in the C# programming language. The mesh used in all the brain examples consisted of 6534 nodes. Visualization of results was implemented with the Gmsh format Geuzaine and Remacle (2009). Table 3 contains the model values. In addition, for the I ion we used the Fitz Hugh Nagumo cell model with rates and constants from the cell ML repository Lloyd et al. (2008). Code available at https://github.com/steppenwolf0/brainSim.

Table 3
Variables values for the simulations.