Background: Ischemic stroke is a major cause of disability and death. Understanding microscale pathophysiology within neurovascular units (NVUs) is critical for developing therapies, yet traditional methods struggle with spatiotemporal heterogeneity. In silico modeling offers a solution for multiscale analysis.

Objective: To develop a 3D in silico model of ischemic stroke integrating NVU morphology, Monte Carlo dynamics, cellular automata, and pharmacodynamics to predict lesion evolution and assess treatments.

Methods: We created a hybrid 3D model featuring: anatomically detailed NVUs; Monte Carlo simulation of oxygen diffusion; cellular automata for cell viability/death rules; a pharmacodynamic module for drug effects. The model was calibrated with in vitro and in vivo (rat) data.

Results: The model replicates: penumbra and infarct core formation (matching experiments); hypoxia and energy deficit gradients in NVUs; neuroprotective drug effects by timing/dose. Sensitivity analysis identified key parameters driving penumbra to necrosis transition.

Conclusion: This 3D in silico platform enables: studying ischemic mechanisms at the NVU level; predicting pharmacotherapy efficacy; optimizing preclinical experiments; studying ischemic injury mechanisms at the NVU level, accounting for stochasticity and discreteness of processes; predicting the efficacy of pharmacotherapeutic interventions in silico; optimizing preclinical experiment design.

The model can serve as a tool for personalizing stroke therapy by simulating individual perfusion parameters and pharmacokinetics.

Ischemia of a specific brain region, resulting from insufficient blood supply, represents one of the most significant and common causes of stroke. It represents a major clinical challenge as stroke can have far-reaching consequences for human health.

Contemporary neurological research demonstrates that ischemia leads to a critical increase in intracellular concentrations of calcium and sodium ions, accompanied by a marked a substantial intracellular acidification occurs. These alterations in the cellular chemical environment adversely affect intracellular energy metabolism, inevitably resulting in cellular damage and, ultimately, cell death.

In the context of developing effective diagnostic and therapeutic approaches for cerebral ischemia, mathematical modeling—particularly the Monte Carlo method-opens new perspectives. This method enables highly accurate simmulate and analysis of the pathogenetic mechanisms underlying ischemia, thereby significantly contributing to a deeper understanding of the processes occurring in the brain during hypoperfusion. It may also serve as a key to developing innovative strategies for stroke prevention and treatment. Ischemic stroke is one of the most prevalent and, unfortunately, extremely severe medical conditions, ranking among the leading causes of disability and high mortality across diverse population groups. This pathological process arises from a critical disruption of normal cerebral blood flow, leading to neuronal cell death, subsequent tissue damage, and consequently, severe impairments in brain function that profoundly affect a person’s entire life.

The aim of this study is to model the pathogenesis of cerebral ischemia in a three-dimensional (3D) brain region using the Monte Carlo method. To achieve this objective, it is necessary to examine the pathogenesis of ischemia, investigate the Monte Carlo method and its applications in medicine, and perform computational modeling followed by result analysis. The work is directed toward improving the understanding of ischemic mechanisms and identifying effective diagnostic and therapeutic strategies. 

Approaches to Modeling Ischemic Stroke

Ischemic stroke is an extremely dangerous medical condition posing serious health risks. It ranks among the leading global causes of death and disability. It occurs due to impaired cerebral circulation, resulting in oxygen deprivation and subsequent neuronal death, which in turn leads to various functional impairments in brain activity.

Among the leading causes of global mortality and disability, ischemic stroke (also known as cerebral ischemia) plays a pivotal role, representing the second-highest cause of death worldwide after cardiomyopathy and affecting tens of millions worldwide. The etiopathogenetic basis of cerebrovascular disorders is cerebral ischemia (~80%) or intracranial hemorrhage (~20%). Ischemia most commonly occurs in the territory supplied by the middle cerebral arteries [1].

In addition to acute neurological deficits directly attributable to the affected brain regions, stroke may also lead to long-term cognitive and psychiatric consequences that are difficult to link specifically to the damaged area. Cognitive decline is a serious issue frequently observed, particularly after lacunar stroke(s), and likely overlaps with underlying small vessel disease [2]. Post-stroke depression is another major concern, developing in approximately one-third of patients with chronic stroke [3]. The scientific community continues to face challenges in understanding the biological underpinnings of certain post-stroke conditions, partly due to a lack of adequate and representative animal models. This significantly hinders research into and development of effective treatments. Moreover, existing animal models often fail to adequately address serious and common complications of chronic stroke, such as obstructive sleep apnea - a condition characterized by periodic cessation of breathing during sleep. This complication is largely overlooked in experimental research, creating gaps in our understanding of the full spectrum of chronic stroke consequences and consequently limiting the development of comprehensive treatment and prevention strategies [4].

Stroke, as one of the most severe and acute conditions in medical practice, affects not only neurological and neuropsychiatric health but also induces significant alterations in immune system function. This phenomenon, known as post-stroke immunodepression, manifests rapidly after stroke onset. Statistical data indicate that approximately 60% of stroke patients experience fever within the first three days following the event. This temperature elevation is a sign of the body’s attempt to respond to the stress induced by stroke and mobilize its defense mechanisms—albeit with delayed and weakened efficacy due to immunodepression [5].

Developing models of ischemic stroke is crucial for understanding the complex interactions within the neurovascular unit, facilitating pathogenesis research, and advancing effective therapeutic strategies. In vitro models, particularly three-dimensional ones, accurately mimic in vivo conditions, enabling precise experimental investigations [6].

As medical science rapidly advances, considerable attention is being devoted to the development and refinement of mathematical models, which play a pivotal role in elucidating the multi-level and complex mechanisms of ischemic stroke. These models are not merely useful—they are becoming indispensable tools for accurately predicting clinical outcomes, thereby significantly enhancing, and optimizing treatment strategies for this serious disease [7].

Mathematical models enable the integration of biological data encompassing numerous parameters and aspects, allowing the creation of extensive, multifunctional databases for analysis. This provides researchers and clinicians with a unique opportunity to simulate diverse clinical scenarios—an invaluable contribution to evaluating the efficacy of various medical interventions and therapeutic approaches.

Mathematical models help unravel the multifaceted biological processes associated with ischemic stroke, ranging from cellular interactions to systemic responses. They facilitate the identification of novel regulatory principles governing disease progression, thereby aiding in the development of targeted therapies [8].

Such models can predict individual patient treatment outcomes based on clinical data, improving prognostic accuracy, and informing clinical decision-making. For instance, a model developed using transcranial magnetic stimulation data effectively stratified patients into groups with favorable versus unfavorable recovery trajectories [9].

Modeling ischemic stroke using the Monte Carlo method represents a modern approach to analyzing and predicting clinical outcomes. This statistically based method allows for highly accurate incorporation of random factors and uncertainties inherent in the development and treatment of ischemic stroke. The Monte Carlo method enables researchers and clinicians to explore diverse clinical scenarios in depth and assess their probabilistic characteristics—critical for designing effective treatment and rehabilitation strategies for patients with this complex condition. This approach not only enhances understanding of disease mechanisms but also significantly improves the quality of care for ischemic stroke patients.

The Monte Carlo method is a unique numerical technique that actively employs random sampling to rigorously model diverse processes and accurately estimate the probabilities of various potential outcomes. The method comprises several key stages, each essential to achieving the goal of analysis and prediction [10].

The first stage involves constructing a mathematical model, which forms the foundational basis of the entire process. At this stage, equations linking input and output variables must be defined, establishing the framework for subsequent computations.

The next step is defining input values. It is crucial to appropriately select probability distributions that adequately represent the stochastic nature of input parameters. This is critical to ensure simulation results closely reflect real-world conditions.

The third stage entails generating a dataset. This requires producing many random samples based on the previously selected distributions, ensuring a sufficiently robust sample size for reliable subsequent analysis.

Finally, the last stage is resulting analysis. Here, the generated data are thoroughly processed to estimate the probabilities of various outcomes, enabling evidence-based conclusions and informed decision-making.

Overall, the Monte Carlo method is a powerful tool widely applied across diverse fields -from financial analysis to engineering-due to its capacity to model complex systems using random sampling. The primary goal of applying the Monte Carlo method in ischemia modeling is to evaluate prognostic factors and disease outcomes [11]. Research shows that this approach enables more accurate prediction of adverse event probabilities based on various clinical parameters, such as fluid-electrolyte balance and patient severity status.

Three-dimensional (3D) mathematical modeling of ischemic stroke offers several key advantages that enhance understanding of pathophysiological processes and improve clinical outcomes. The main aspects underscoring the importance of this approach include:

1. Detailed Representation of Anatomy and Physiology

Three-dimensional modeling enables more accurate reproduction of the brain’s complex anatomical structures, including the vascular network and neuronal connectivity. This is essential for understanding how ischemia affects different brain regions and how it propagates. 

2. Dynamic Process Modeling

Three-dimensional (3D) models can account for dynamic changes in cellular and tissue states over time. This enables the investigation of processes such as neuronal depolarization, shifts in ion concentrations, and metabolic reactions, which are impossible to replicate in static models.

3. Predicting Clinical Outcomes

The application of 3D models for predicting outcomes in ischemic stroke allows for the integration of numerous factors, including the degree of tissue damage, the rate of reperfusion, and the effects of therapeutic interventions. This integrated approach facilitates a more accurate prognosis for the likelihood of functional recovery and patient survival.

4. Optimization of Therapeutic Strategies

Modeling can assist in the development of personalized treatment approaches by allowing the evaluation of the effectiveness of various therapies depending on the specifics of each case. This is particularly relevant in the acute phase of a stroke, where time is of the essence.

5. Enhancement of Educational Processes

Three-dimensional models can be used for training medical professionals, providing them with the opportunity to visualize the complex processes occurring in the brain during an ischemic stroke. This improves specialists' understanding and preparedness for clinical practice.

6. Investigation of Novel Treatment Methods

3D-level modeling opens up possibilities for testing new treatment methods and technologies, such as stenting or thrombolysis, in controlled conditions before their practical application.

Mathematical modeling of ischemic stroke at the three-dimensional level represents a crucial tool for enhancing the diagnosis, treatment, and rehabilitation of patients with this condition. It enables a deeper understanding of pathological mechanisms, optimization of clinical approaches, and development of novel therapeutic methods, ultimately contributing to reduced stroke-related morbidity and mortality [12].

In this scientific work, we rely on a number of fundamental principles that form the basis of the model we employ [15-24]. These principles represent key aspects that define the structure and functionality of the system under investigation. Specifically, we focus on the discrete states of the morphofunctional units of the brain, known as neurovascular units (NVUs). It should be emphasized that within the framework of this model, we do not delve into the details of the interaction mechanisms between these units, leaving them outside the scope of our current research.

Furthermore, it is important to note that our model does not address the process of forming the primary zone of cerebral blood flow deficiency. In the context of our study, this zone is considered as initially given, which simplifies the task and allows us to focus on other aspects of the model.

The main system variables we investigate are related to spatial dynamics and include the specific content of NVUs residing in each of the discrete states. These variables correspond to the total number of NVUs located within a specific, finite-sized area of brain tissue. The study of these variables allows us to better understand the processes occurring within the system under investigation.

In addition to the above, it has been decided not to employ a continuum model, which involves averaging of cell characteristics and can introduce errors in describing intercellular relationships. Instead, a discrete model is applied, where the neurovascular unit (NVU) serves as the minimal functional unit. The research utilizes principles for constructing such a model that are highly efficient and do not require significant computational time. This makes the model not only theoretically sound but also practically applicable for further research in this field.

This work focuses on studying a three-dimensional model that provides a simplified, yet highly informative representation of the complex processes occurring in cell culture. The research concentrates on the behavior of a specific section of cellular tissue, whose boundaries are determined by our choice and defined solely by the area of interest. It is important to emphasize that the size of the studied fragment is negligible compared to the total volume of the brain, allowing us to confidently consider this area as part of an infinite two-dimensional space where Euclidean metric is applied.

Particular attention in the study is given to the propagation velocity of various substances within neurons and astrocytes, as well as the transmission speed of the membrane potential along cell membranes. These processes occur at a remarkably high rate, significantly exceeding the speed of ionic transport across the cell membrane. Additionally, it is important to note that within the area supplied by a single arteriole, we observe homogeneity in parameters such as partial pressure, concentration of various biologically active substances (including ions and oxygen), and enzyme activity. Based on this, it is hypothesized that the space within the neurovascular unit (NVU) can be considered homogeneous, which constitutes a key aspect of our hypothesis.

Finally, at the initial observation time point, denoted as t = t0, it is assumed that a specific, albeit small, proportion of cells had been subjected to acute ischemia. This assumption plays an important role in the subsequent analysis and interpretation of the obtained data.

2.1. Description of the Evolution of Lattice Fragment States

Under certain conditions, a transition between different states of the system becomes possible. For the purpose of the model, a lattice fragment consisting of N × N elements is considered, with appropriate boundary conditions applied. Within the Markov approximation, the evolution of the probability distribution of the fragment states is described by the master equation (ME):

t - time (s); S - the state of the lattice fragment at time t, defined by the occupation numbers si of all nodes, S = (si), i=1,2,…,N2, si ∈ {1,2,3,4};  PS(t) - the unconditional probability of the system being in state S at time t; λ(S′ → S, t) is the transition rate from state S′ to S at time t, with initial conditions Ps(0) = ????????0

The intensity of transition from one state to another depends on the total rates of all possible elementary processes occurring within the considered lattice fragment and is directly influenced by its’ previous configuration. In this article both single-site and two-site processes occurring between nearest-neighbor nodes are examined. Their rates are calculated according to transition-state theory as follows:

α ∈ {1; 1−; 4} and β ∈ {2; 3; 5; 6; 7; 8} - the numbers of elementary processes (transition stages); i, j ∈ {1, 2, …, N²}  - the numbers of lattice node; δij,β(t,a,b) = P(si = a, sj = b) + P(si = b, sj = a/a ≠ b); a and b correspond to the interacting stages of the elementary events with indices α (β); kα, kβ –constants of rates of respective transition stages.

Because the master kinetic equation (1) constitutes a system of high dimensionality, approximate numerical methods are typically employed for its solution [13]. Individual trajectories of the system in the state space can be generated using the Monte Carlo method. In this study, a dynamic Monte Carlo algorithm with a variable time step was used to simulate the processes occurring on a planar brain region represented by a lattice fragment comprising N × N centers. The algorithm proceeds through the following steps:

a) It is assumed that at the initial moment t = 0, the system is in state S;

b) The rates of all possible surface processes and the total rate of state change are computed:

c) the moment when the system exits the current state is determined, ∆t, is determined as a random variable:

E - a random variable uniformly distributed over the interval (0…1);

d) One of the surface processes available in the current state is selected with a probability proportional to its rate, the state of the lattice fragment is updated accordingly, and the procedure returns to step (b).

Such a computational framework is referred to as a microscopic stochastic or simulation model. In this simulation model, the dynamics of the selected lattice fragment are described in terms of the occupation numbers of adsorption centers, thereby representing processes at the level of neurovascular units (NVUs) in brain tissue.

2.2. Elementary Processes

To designate the states of individual cells, letter-based mnemonic codes were used:

   H — the healthy state of a neurovascular unit (NVU);

   D — the dead state of the NVU (resulting from necrosis or apoptosis);

   I — the ischemic state of the NVU;

   F — the state of pharmacologically active substance (PAS) accumulation within the NVU.

   F — the state of accumulation of a pharmacologically active substance (PAS) in NVU.

Figure 1: Rules of the cellular automaton model (a); in the presence of a pharmacologically active substance (b)

In the model, the state of the lattice is updated according to a set of rules that refers to real physiological processes occurring in both healthy brain regions and ischemic areas. The sequence of elementary surface processes is generated using the dynamic continuous-time Monte Carlo method described above.

	Transition of a cell from the healthy state to the ischemic state
	Transition from the ischemic state to the cell death state
	Transition from the ischemic state to the healthy state
	Transition of an ischemic cell to the state of PAS accumulation
	Transition from the PAS accumulation state to the healthy state
	Return from the PAS accumulation state to the ischemic state

Table 1: Single-site processes of the model.

	Transition from the ischemic state to the healthy state in the presence of a neighboring healthy cell
	Transition of a healthy cell to the ischemic state in the presence of an ischemic cell
	Transition of a healthy cell to the ischemic state when adjacent to a dead cell
	Transition of an ischemic cell to the dead state when adjacent to a dead cell
	Transition of an ischemic cell to the dead state when adjacent to another ischemic cell

Table 2: Two-site processes of the model

2.3 The relationship between the developed model and the biological model.

As previously mentioned, the unit of the cellular automaton in the developed model corresponds to the state of a neurovascular unit in the brain. Ischemic stroke is associated with damage to small, penetrating cerebral arteries measuring 40 to 200 microns in diameter, which supply blood to deep regions of the brain [14]. NVU is a structural and functional unit comprising vascular and cellular components in the vicinity of the terminal branches of these penetrating arteries, including arterioles and capillaries.

According to current research, the occlusion of a single small artery results in the so-called lacunar infarction, which corresponds to a region of brain tissue supplied by that artery. Lacunar infarcts manifest as small areas of tissue destruction measuring between 3 and 15 millimeters in diameter, or as small cavities that may contain multiple smaller foci. For the most part, these are small infarcts that correspond to the area of blood supply of an artery. Additionally, healed or absorbed microbleeds are common and can only be detected through modern, high-precision research techniques and morphological examination. Using radiation diagnostics, lacunar infarcts have corresponding cerebrospinal fluid.

At the same time, other researchers define foci up to 1 centimeter in diameter as small strokes, and foci around 2 centimeters in diameter as large heart attacks, with an average of small foci. Therefore, it seems reasonable to consider a 0.1 centimeter stroke focus as the biological scale corresponding to the modelwe are developing to study the model proposed in this thesis. Along a segment equal to this diameter, approximately 33 high-resolution elements can be placed. Consequently, the 30×30 grid size of the developed model can be considered consistent with the biological characteristics of the modeled object.

The program for modeling ischemic stroke was implemented in Python using standard libraries, including numpy, matplotlib.pyplot, and scipy.integrate.

The numpy library's linspace function was used to create a sequence of evenly spaced data points on a numeric axis within a specified interval. The same library's meshgrid function was used to generate a two-dimensional grid of x-y coordinates based on specified numerical ranges. This allows for the transformation of one-dimensional arrays into two-dimensional matrices, providing a universal grid of x-y coordinate pairs required for various functions.

Python's odeint function, available through the scipy.integrate library, was utilized to solve ordinary differential equations. The odeint() function returns an array with length equal to the number of time points (len(t)) multiplied by the initial condition vector length (len(y0)). The odeint () function offers numerous options that allow for customization of its operation. Differential equations described in Chapter 2 were solved utilizing this function.

Pyplot is a Python module within the Matplotlib data visualization package. It automates the creation of axes, shapes, and other elements, eliminating the need for manual intervention. The module provides a set of command-style functions, allowing users to interact with Matplotlib in a manner similar to MATLAB. All graphs generated in this study were created using the Pyplot module.

As a result of testing the program with the assignment of empirically selected coefficients, data were obtained on the evolution of the states of the lattice of the cellular automaton of the model for the spread of stroke in a region of the brain tissue. The results of one run are presented in Fig. 3 as a graph of the time distribution function of the specific number of each cell type in the stroke focus.

Figure 3: The result of testing the program with the appointment of empirically selected coefficients. On the vertical axis - the specific number of cells in each state.

Figure 4: 2D Visualization of the initial state of the stroke focus by the simulation program. Dead cells are black, ischemic cells are red, healthy tissue cells are blue.

Figure 5: 2D Visualization of the infarction zone spread and penumbra formation by the simulation program. Dead cells are black, ischemic cells are red, healthy tissue cells are blue. Penumbra is the zone between blue and black density.

Here is the initial focus of ischemia with a lesion of 10% of the final size of the focus. In the future (Fig.5), as a result of modeling the interaction of neurovascular units, the infarction zone spreads and the formation of the border zone of the penumbra, experiencing, on the one hand, the influence of necrotic and ischemic focus cells, and on the other, relatively intact tissue.

As a result of the program, a three-dimensional (3D) visualization of the intensity of ischemic stroke was obtained, Figure 6 (6.1-6.3).

Figure 6.1: 3D Visualization Spatial model of ischemic stroke intensity. Initial state.

Figure 6.2: 3D Visualization Spatial model of ischemic stroke intensity. The intermediate state.

Figure 6.3: 3D Visualization Spatial model of ischemic stroke intensity. The formed focus of stroke.

Within the framework of this project, a model has been implemented based on the application of a set of ordinary and partial differential equations. This is a crucial aspect in accurately reproducing complex biological processes.

The aim of this work was to develop software that can simulate the process of ischemic stroke. In order to do so, we used the parameters that constitute the reference set of values. As a result, we obtained a three-dimensional model of ischemic stroke.

An innovative computational platform has been developed. We have created the first 3D in silico model of ischemic stroke that comprehensively integrates: anatomical organization of neurovascular units (NVUs); stochastic dynamics of substance transport (Monte Carlo method); discrete cellular processes (cellular automata); pharmacodynamic interactions.

The model reproduces key pathophysiological phenomena. Computational experiments confirmed that the model adequately simulates: formation of the penumbra and infarct core with realistic spatiotemporal dynamics; gradients of hypoxia and energy deficit within NVUs; dose  and time dependent effects of neuroprotective drugs.

Critical parameters of pathogenesis have been identified. Sensitivity analysis revealed key factors determining the transition of cells from penumbra to necrosis. This opens opportunities for targeted intervention in the «bottleneck» stages of the ischemic cascade.

The model has multidisciplinary value. The platform serves as a tool for: fundamental research into microscale mechanisms of ischemic injury; preclinical assessment of pharmacotherapeutic strategies in silico; optimization of biological experiment design through computational hypothesis screening.

Personalized medicine prospects are opened. Thanks to the ability to calibrate the model to individual perfusion and pharmacokinetic parameters, it can become a foundation for personalized stroke therapy — predicting treatment response for a specific patient prior to administration.

The potential of hybrid approaches is demonstrated. The successful integration of Monte Carlo methods, cellular automata, and pharmacodynamics confirms the effectiveness of hybrid computational strategies for modeling complex biomedical systems where stochastic, discrete, and continuous processes coexist.

Future development directions are defined. Promising tasks include: validation of the model using clinical data; expansion of the pharmacodynamic module to simulate drug combinations; implementation of machine learning for automated parameter calibration.

In summary, the proposed 3D in silico model represents a scalable and adaptive tool for studying ischemic stroke and developing novel therapeutic strategies.

The work was supported by the grant of the Ministry of Science and Higher Education of the Russian Federation No. 075-15-2020-792 (Unique identifier RF----190220X0031).

Authors' contributions: 

Dr. Valentin V. Fursov - idea, research management, mathematical modeling per se, text, editing supervising; MSc Alexander V. Ananiev – 2D mathematical modeling, program code, STDT. Ksenia A. Popova – 3D modeling program code, Dr. Ilia V. Fursov – biomedical aspects, MSc. Ivan V. Fursov – neuroscience aspects, STDT. Anikey V. Fursov – modeling, coding & paper preparing assistance.

Acknowledgments: 

Prof. Dr. Dmitry A. Kuznetsov, Prof. Dr. Eleonora M. Koltsova, Dr. Ivan I. Mitrichev, Moscow, Russian Federation, are to be sincerely thanked for their kind assistance with preparation of a manuscript, consultations and for their encouraging comments on this work.

Conflicts of interest 

No conflicts of interests of any sort beyond.

Informed Consent: none to declare

Data Availability:

None to declare