In Silico Studies on Pharmacokinetics and Neuroprotective Potential of ²⁵Mg²⁺: Releasing Nanocationites - Background and Perspectives

Abstract

Sharp blood circulation disorders are known for their capability to promote such abundant and hardly treatable pathologies as myocardium infarction and the ischemic brain stroke (“insult”). Noteworthy, the stroke — related brain tissue metabolic damages involve an essential ATP deplete clash along with a suppression of brain specific nucleotide — associated kinases and ATP synthase, both Mg2+ — dependent complex enzyme “machineries”. This itself makes the latter’s a legitimate target for some advanced pharmaceuticals as long as the drug — induced overstimulation of corresponding enzymatic activity is the case. Thus, magnetic isotope effects (MIE) of the nuclear spin possessing paramagnetic 25Mg2+ ions might modulate the brain creatine kinase, alfa-glycerophosphate kinase and pyruvate kinase catalytic activities in a way of a remarkable ATP hyperproduction required to compensate the hypoxia caused acute metabolic breakdown. To realize the Magnesium-25 pharmacological potential, a low-toxic amphiphilic cationite nanoparticles were introduced lately. Particularly, the Magnesium — releasing porphyrin-fullerene nanoadduct (cyclohexyl-C60-porphyrin, PMC16) has been proposed to meet expectations dealing with a targeted delivery of 25Mg2+ towards the brain ischemia surrounding areas. In order to optimize a multi-step [25Mg2+]4PMC16 preclinical trial scenario, the In Silico algorithms are to be developed and analyzed. In this study, these algorithms are in a focus with a special emphasize on a novel combination of slightly modified Gompertzian equation systems and a non-Markov population dynamics concept. This In Silico approach takes into account some literature-available patterns of brain hypoxia pathogenesis, the resulted simulation model could be considered as a promising tool for further research on experimental nanopharmacology of the ischemic stroke.

1. Introduction. Formulation of the problem

Research in the field of development of promising drugs for the treatment of ischemic stroke based on nanoparticle carriers of cations of paramagnetic isotopes of divalent metals is at an early stage. Such studies are compounded by the lack of descriptions of relevant mathematical models of ischemic stroke in scientific literature, as well as drug-specific models of the pharmacokinetics of targeted delivery of PMC16 nanoparticles. All this is required for the development of In Silico instruments for preclinical and clinical studies of the neuroprotective potential of [²⁵Mg²⁺]₄PMC16 in the treatment of ischemic stroke [1, 2, 3, 4] (Figure 1).

Figure 1.

Prospects for the use in neurobiology of “smart” nano-cationites based on porphyrin adducts of fullerene C60.

Creation of In Silico - algorithms for optimization of multistage scenarios of preclinical trials of [²⁵Mg²⁺]₄PMC16 in experimental nanopharmacology of ischemic stroke represents a completely new, complex, innovative and challenging interdisciplinary problem.

Speaking of the direct and clear practical benefits which are supposed to be gained from the appropriate use of mathematical modeling in specifying the plan of preclinical anti-hypoxia medicines research, this requirement is undoubtedly essential for optimizing this plan. Notably, an applied pharmacological potential of such a peculiar In Silico simulation approach might be taken as a “hopeful pullout” for coming up with a novel element in a preclinical trial strategy for prevention of metabolic breakdown in brain ischemia and/or correction based on the administration of paramagnetic bivalent metal isotopes released and delivered by amphiphilic nano-cationites belonging to the superfamily of PMC16 (C60-porphyrin) nanoparticles [5].

Noteworthy, a so-called “sovereign trend” in computational modeling of pharmacological processes within the current preclinical trial paradigm has already made a significant impact on preclinical trial design in experimental neurology and neuropharmacology [6, 7, 8]. This correlates with the PubMed statistics showing a remarkable increase in the number of publications on the above-specified issue [6].

2. Methods

Simulation of processes occurring in biological objects and systems is necessary to optimize algorithms for preclinical and clinical studies in pharmacology. These tasks are solved using in vitro, in vivo, In Silico models. Modeling is one of the leading research methods of this kind. The variety of processes in a living organism is so great that it is almost impossible to get a detailed and complete understanding of the behavior of a living system. In view of this, the development of new treatment methods, diagnosis, pharmacy, etc., requires the modeling of objects of appropriate research. Any type of modeling consists of replacing the investigated object (process, phenomenon) with a model, which is a semblance of a real object (process, phenomenon). At the same time, such an object representing the model is consciously perceived as simplified. However, it is vital that it retains the main, most essential properties for research, which are available for a real object (system, process, phenomenon).

Modeling is a method in which the study of its model replaces the study of a complex object (process, phenomenon). Accordingly, such an object (process, phenomenon) itself, which resembles the real object, but has been deliberately simplified, is called a model.

Any scientific research method, including both theoretical and empirical, is based on the idea of modeling.

In this work, we will adhere to the classical algorithm for constructing mathematical models adopted in biophysics.

The main stages of modeling can be summarized as follows:

1. Primary collection of information about the object of modeling: about its properties, processes occurring in it, patterns of behavior under various external conditions.

2. Formulation of the problem. The goal of the study and its main tasks are formulated. It is determined what new knowledge should be obtained after the research has been conducted.

3. Substantiation of basic assumptions. It is necessary to determine the characteristics of the object that are insignificant for solving the research problem, which can be neglected.

4. Creating a model, researching it.

5. Checking the relevance of the model to the object under study.

3. General task structure

The difficulty lies not only in the fact that in the domestic and foreign literature, there are no relevant mathematical models of ischemic stroke, and they need to be created almost from scratch, but also within the problem itself, which arises from the necessity not only to develop but also to align mathematical models of several processes mutually:

1. The process of necrosis of brain tissue as a result of ischemia and related phenomena (apoptosis, toxicosis, edema, etc.) in the absence of pharmacotherapy;

2. Pharmacokinetics (i.e., delivery) of the [²⁵Mg²⁺]₄PMC16 drug to the desired area of the brain and distribution throughout the tissue;

3. Model of the process of the effect of the drug [²⁵Mg²⁺]₄PMC16 on the synthesis of ATP and the prolongation of the life cycle of cells, which are subjected to hypoxic conditions, but have not lost their viability;

4. The recovery process (reperfusion, [neuro]glialisation, regeneration) of the functions of living ischemic cells as a consequence of the pharmacotherapy of ischemic stroke with the drug [²⁵Mg²⁺]₄PMC16.

5. Phagocytosis.

6. Other processes.

The general structure of the problem is shown in (Figure 2).

Figure 2.

General structure of the task of developing In Silico pharmacokinetics and neuroprotective potential of fullerene-porphyrin nano-cationites carrying ²⁵Mg²⁺ in relation to the pathogenesis of ischemic stroke.

At the same time, the In Silico development process is implemented in several stages:

1. Formation of a hypothesis and a general structural model In Silico.

2. Acceptance of initial constraints and simplifications.

3. Formation of hypotheses and primary models of processes.

4. Mathematical modeling of individual processes.

5. Combining mathematical models of individual processes into a system of differential equations.

6. Search for optimal solutions to the system of differential equations (In Silico level I).

7. In Silico level II development: algorithmization, IT programming, debugging, testing.

8. Use of level II In Silico for prediction and simulation of in vivo.

9. Raising accuracy, development, improvement.

10. Development of “smart” In Silico level III based on artificial intelligence neural networks: Design of neural network architecture, debugging, training, testing.

11. Use of “smart” In Silico level III as a predictive tool for preclinical trials.

12. Laying the foundations of smart In Silico level IV architecture as a predictive tool for clinical trials.

4. Mathematical modeling of ischemic stroke. Hypothesis formation

Mathematical modeling of ischemic stroke is a complex task in itself. Several pathogenetic subtypes (atherothrombotic, cardioembolic, lacunar, hemodynamic and microcirculatory) have been highlighted. Accordingly, the mechanisms of occurrence and development of the disease also differ. All this significantly complicates the modeling of the development of this disease at the level of hypotheses and primary algorithms laid down in the In Silico process. This nosology of ischemic stroke complicates approaches to the formation of primary algorithms for the process of occurrence and course of ischemic stroke, as well as solving problems of mathematical formalization.

It is customary to distinguish 4 stages [9] of ischemic stroke (Figure 3):

1. Terminated (3–5 days)

2. Most acute (7–10 days)

3. Acute (up to 1 month)

4. Early recovery (up to 6 months)

5. Late recovery (from 6 months to 1 year)

6. Long-term (over 1 year)

Figure 3.

Graphic model of ischemic stroke, the development of ischemic penumbra and the spread of the damaged brain tissue due to ischemia and blockage of the vascular bed branch.

Of these, the acute phase is characterized by the most severe course and high mortality, which is why the study of patterns of the disease in the acute phase (first five days of the disease) is the most significant. During this period, the right therapeutic tactics can bring the maximum result. During this period, the drug [²⁵Mg²⁺]₄PMC16 should have the maximum positive effect on the dynamics of the course of the disease.

In the development process in silico, we will build on the simplest concepts of the course of the disease, introducing some assumptions and simplifications. According to the need to improve the model, these simplifications will subsequently be removed or replaced by more complex designs. Then the problem of merging the models of the various processes as mentioned earlier into a single structure needs to be solved.

Based on this, we will simplify the growth of an ischemic stroke as a process of necrotic death of cells of normal brain tissue caused by a stop in the blood supply to its part due to blockage of a part of the vascular bed that feeds the brain. Consequently, hypoxia begins to develop in the area of the brain, the blood supply of which was carried out from the clogged vascular network. As a result of hypoxia, cells that lack oxygen begin to die off gradually. A process of necrosis of this part of the brain’s tissues takes place (Figure 3).

The process of necrotic cell death is stretched out in time. For this reason, there is a transition region between the infarct nucleus and healthy tissue (ischemic penumbra), in which the functions of the cells are disrupted, but they remain viable. It is vital to note that penumbra cells’ death is reversible and progresses more slowly than in the infarct nucleus, within a few hours [10]. This “therapeutic window” (at least 3–6 hours) provides time for diagnosis and treatment measures aimed at restoring nerve cells, limiting the area of damage and reducing neurological consequences.

Thus, in order to solve the general problem: the in silico development of the neuroprotective potential of [²⁵Mg²⁺]₄PMC16 nanoparticles, it is necessary to develop in silico of ischemic stroke and combine this model with the pharmacokinetics of [²⁵Mg²⁺]₄PMC16 as a drug. At the same time, it is necessary to consider the particularities of targeting (targeted delivery) of [²⁵Mg²⁺]₄PMC16 nanoparticles in the infarction zone, more precisely in the area of ischemic penumbra.

5. Acceptance of initial restrictions and simplifications

Starting the mathematical modeling of ischemic stroke at the cellular level, we will accept a number of assumptions and simplifications, namely:

1. We consider the brain to be homogeneous in composition, structure and cell type. Cell differentiation is neglected;

2. We do not consider any other causes of cell death other than necrosis as a result of hypoxia;

3. We do not consider the processes of removing the decay products of dead cells and toxins that occur in the process of brain activity;

4. We do not take into account the consequences of brain tissue necrosis, such as edema, an increase in the volume of necrotizing tissue, leading, in particular, to a spike in intracranial pressure, deformation of brain structure, etc.

5. We will assume that brain volume and the number of its cells before, after and during the development of a stroke remain unchanged;

6. We will assume that there is a certain “point of no return” for brain cells in the process of hypoxia, after which the process of cell restoration (regeneration) is not possible. This factor affects the size of the ischemic penumbra zone;

7. We will assume that during the period considered in silico, the general blood supply to the brain and the body as a whole does not stop (the patient does not die);

8. We will assume that no outside interference (both therapeutic and surgical) is carried out in the body during the period under consideration.

9. In this way, the general model of the dynamics of the development of ischemic stroke in time will look as follows, presented below (Figure 4):

Figure 4.

Time model of progression of ischemic stroke. In the figure: t0 is the moment when blood supply stops and hypoxia begins, t1 is the onset of cell necrosis, t2 is the start of recovery processes, t2 is the moment of cessation of brain cell death, t3 is the moment when hypoxia ends (blood supply is restored), t4 is the end of cell necrosis, t5 is the moment of completion of the post-stroke recovery process.

6. Formation of hypotheses and primary models of processes

Then, in light of the assumptions made, for the purposes of mathematical modeling, the hypothesis of the scenario of ischemic stroke progression will look in the following manner.

Until time $t_0$ (pre-stroke phase), the brain functions normally, the blood supply is not disrupted. At time $t_0$ instantaneous (additional) cessation of blood supply occurs, and the process of necrotization of the brain tissue begins. We believe cell death in hypoxia is also stretched over time and occurs in stages. That is, up to a specific time period $t_1$, the cell is still capable of recovery, and we can talk about the existence of a “point of no return”, after which cell recovery is impossible. We will take this into account later. At time $t_2$ recovery processes in the body are activated, due to the entry of blood supply from the brain area that has not undergone necrosis. Thus, we can talk about the presence of two “counter” processes in the model of ischemic stroke: the process of dying (necrotization) of cells and the process of restoring the functions of brain cells (Figure 5).

Figure 5.

Counter processes and the resulting process of ischemic stroke growth.

Based on this, we can distinguish the following phases of ischemic stroke:

1. Pre-stroke phase, t≤ $t_0$,

2. The phase of brain cell death (necrotic phase), $t_0$ ≤ t ≤ $t_4$;

3. Brain cell recovery phase, $t_2$ ≤ t ≤ $t_3$;

4. Post-stroke phase, t ≥ $t_5$.

7. Mathematical modeling of individual processes

Now, we can write in general algebraic form the first simple equation, which we will call the “hypoxia equation”, describing the amount of GM tissue that has undergone hypoxia at time $∆t=t-t_0$, taking $t_0=0$ as the starting point of the hypoxia period. ($∆t=t$):

$$Q_g=v_g\ast t$$E1

Where $Q_g$ – the amount of GM tissue that has undergone hypoxia,

$v_g$ – specific rate of hypoxia (death) of the brain tissue,

t – time elapsed since the end of the blood supply to the brain area.

In this case, the amount of necrotizing tissue can be represented as its volume V, mass m, or the number of dead $N_d$ cells. At the same time, the structure of the formula will not change. Analogously to formula Eq. (1), we write down the necrosis formula for the number of dead $N_d$ cells:

$$\begin{array}{c}{N}_d\left(t\right)=\delta \ast t\end{array}$$E2

Where $N_d\left(t\right)$ – the number of dead cells since the end of the blood supply,

$\delta$ – specific rate of cell necrosis (death),

t – the time elapsed since the end of the blood supply.

Or, in differential notation:

$$\begin{array}{c}\frac{d{N}_d\left(t\right)}{\mathit{dt}}=\delta N_d\end{array}$$E3

Similarly, let us write in differential form the formula for the regeneration of brain cells in the ischemic zone.

$$\begin{array}{c}\frac{d{N}_R\left(t\right)}{\mathit{dt}}=\alpha N_R\end{array}$$E4

Where $N_R\left(t\right)$ – number of recovered cells,

$\alpha$ – the specific rate of cell repair.

Under our assumptions, the conservation law (balance formula) will look like this:

$$\begin{array}{c}{N}_{\mathit{tot}}=N_n\left(t\right)+N_g\left(t\right)\end{array}$$E5

Where $N_{\mathit{tot}}$ – the total number of brain cells,

$N_n\left(t\right)$ – the number of living normal cells not exposed to hypoxia,

$N_g$(t) – the number of hypoxidated cells and those at different stages of necrosis (dying).

Assuming that the volume of the brain does not change, and due to the assumption of its homogeneity and the absence of processes for removing dead and damaged cells, the number of brain cells also does not change, we can write:

$$\begin{array}{c}{N}_{\mathit{tot}}=\mathit{\text{const}}\end{array}$$E6

It is obvious that the number of hypoxidated cells consists of dead cells that are no longer subject to restoration and regenerated cells, since the processes of dying and restoration in the population of hypoxidated cells proceed simultaneously.

$$\begin{array}{c}{N}_g\left(t\right)=N_d\left(t\right)+N_R\left(t\right)\end{array}$$E7

Then equation Eq. (7) can be rewritten as follows:

$$\begin{array}{c}{N}_n\left(t\right)+N_d\left(t\right)+N_R\left(t\right)=\mathit{\text{const}}\end{array}$$E8

Based on our assumptions, and considering that:

$$\begin{array}{c}{N}_{\mathit{tot}}=\rho \ast m\end{array}$$E9

Where $\rho$ – cell density per unit of brain mass, m is brain mass.

The balance formula Eq. (8) can be rewritten as follows:

$$\begin{array}{c}{N}_n\left(t\right)+N_d\left(t\right)+N_R\left(t\right)=\mathit{\rho m}\end{array}$$E10

where in our assumptions $\rho ,m=\mathit{\text{const}}$.

The derived formula will be useful in the future for clarifying In Silico by in vitro because in biological models and in humans, the volume and mass of the brain, and therefore the number of brain cells, are different. If no distinction is made between healthy normal cells, which have not been exposed to hypoxia, and restored cells, their sum, i.e. the number of functioning cells, can be designated as N(t). Then the balanced equation will take the following form:

$$\begin{array}{c}N\left(t\right)=\mathit{\rho m}-N_d\left(t\right)\end{array}$$E11

Let us consider a spatial model of progression of necrosis and restoration of brain tissue in the process of hypoxia (Figure 6), where we distinguish 3 fundamentally different areas:

Figure 6.

Brain tissue at the border of the hypoxic zone and the normal zone. Area of normal tissue that has not undergone hypoxia (normal area); 2. Area of dead tissue (area of necrosis); 3. Area of “semi-dead” tissue (ischemic penumbra).

Typically, the blood supply to a particular area of the brain is carried out from all sides via the branched vascular network that runs through the brain from the arteries. In case of a local stop of blood supply, the brain will be fed through the active vascular branches located on the periphery of the hypoxic zone. Thus, the cells closest to the focus of the stroke will experience the greatest shortage of blood supply (the focus of the stroke), and the cells located on the periphery will be exposed to the opposite effect (ischemic penumbra).

So, we can imagine that the wave of necrotization, i.e. area 3 (Figure 6), from the moment of local blood supply failure (the beginning of hypoxia), spreads from the focus of the stroke (the place of blockage of the vascular bed) to the periphery. The recovery process goes in the opposite direction: from the periphery to the center of ischemia. After a short period of time after a local stop of blood circulation and the occurrence of hypoxia, necrotization of brain tissue begins. This area of necrotization, expanding, capturing more and more arrays of healthy cells, quickly spreads to the periphery until it meets the area that receives sufficient nutrition from the vascular branches that are not affected by the stroke. Meanwhile, during the development of a stroke, the body starts the recovery processes, increasing the blood supply to the healthy branch. At some point, the “necrotizing wave”, which can be called the “stroke front”, reaches the brain’s area that receives sufficient blood supply from the neighboring unclogged vascular branches. By this time, adaptive processes have already been activated in the body, and through the neighboring non-clogged vascular branches, an increased blood supply is carried out, compared with the norm, sufficient to restore the functions of the brain tissue in the nearby ischemic penumbra. Here, the functions of the cells that are still capable of this restoration are restored. These cells are put back into operation. In this zone, the recovery process begins to dominate over the necrotic one, and the “stroke front” begins to move in the opposite direction (to the focus of the stroke) under the influence of the recovery process. This recovery process is carried out in the ischemic penumbra until the cells capable of recovery are completely exhausted. This process affects all cells for which recovery is still possible. The cells that have finally died are phagocytized, and their remains are removed from the body through the active vascular networks. The released space is filled with connective or adipose tissue.

This complex process can be considered as two opposite “counter” processes: the process of necrosis and the process of recovery. As mentioned earlier, the recovery process is oppositely directed to the process of necrosis and passes with some delay in time.

They can be represented on average as one, a process that is carried out with a certain specific variable total speed, which consists of the specific speeds of necrotization and recovery.

We can write this circumstance as the sum of the corresponding functions as:

$$\begin{array}{c}\epsilon \left(t\right)=\alpha \left(t\right)+\delta \left(t\right)\end{array}$$E12

Then, the differential equation Eq. (3) is transformed into the following form:

$$\begin{array}{c}\frac{d{N}_d\left(t\right)}{\mathit{dt}}=\epsilon N_d\end{array}$$E13

Where $\epsilon \left(t\right)$ – the specific total speed of the ischemic stroke process, is a rather complex function of many variables.

By rewriting equation Eq. (13) in the form:

$$\begin{array}{c}d{N}_d\left(t\right)=\epsilon N_d\mathit{dt}\end{array}$$E14

The solution of this equation in general form can be obtained by integrating over t in the range [t₀, t]:

$$\begin{array}{c}{N}_d\left(t\right)=\underset{t_0}{\overset{t}{\int }}\epsilon N_d\mathit{dt}\end{array}$$E15

At $\epsilon =\mathit{\text{const}},$which can be understood as the average specific rate of development of necrosis in stroke. The trivial solution to this equation for $N=N_{0,}$ is the exponent:

$$\begin{array}{c}{N}_d\left(t\right)=N_0{e}^{\mathit{\epsilon t}}\left(\mathrm{t}\ge 0\right)\end{array}$$E16

Where $N_{0,}$ – the number of brain cells that have undergone hypoxia.

In the case of $\epsilon <0$ formula Eq. (17) represents a specific case of a process where necrosis dominates, and eventually, all brain cells that have undergone hypoxia die. In this case, the formula, which we obtained generally reflects the dynamics of necrosis, which leads to the complete death of the population of brain cells, coinciding with the population dynamics according to Malthus (Figure 7).

Figure 7.

Model graph of brain cell necrosis.

Furthermore, substituting the obtained formula Eq. (16) into the balance equation Eq. (11), we have:

$$\begin{array}{c}N\left(t\right)=\mathit{\rho m}-N_0{e}^{\mathit{\epsilon t}}\left(\mathrm{t}\ge 0\right)\end{array}$$E17

Where $N\left(t\right)$ – the number of living, functioning (including cells that have undergone hypoxia, but have regained their functions) brain cells during the course of an ischemic stroke,

$\rho$ – brain cells density,

m – brain mass,

$N_0$ – the number of brain cells that underwent hypoxia during the [t₀, t] course of ischemic stroke,

$\epsilon$ – the specific average rate of ischemic stroke during the period from the beginning of hypoxia to recovery.

The resultant formula is a level I mathematical model obtained with the maximum simplification of the process of ischemic stroke. For large $\epsilon >0$ it has no biological meaning. For small positive $\epsilon \ge 0$, ε → 0, it exponentially approaches from the value of $N_0$ to the value of $\mathit{\rho m}$, i.e. it mainly describes the recovery process after the acute phase of an ischemic stroke. The model does not describe the initial phase of an ischemic stroke in sufficient detail, resulting from the simplification made when obtaining the formula (17). Additionally, the model does not clearly express the phases of ischemic stroke development.

Accordingly, in the most simplified form, we have obtained a mathematical model of level I, which reflects the essence and the general dynamics of ischemic stroke in borderline conditions, but does not take into account many parameters and, accordingly, does not have the required predictive accuracy. According to the level I model, the graph of the development of ischemic stroke is demonstrated in (Figure 8).

Figure 8.

Graph of the growth of ischemic stroke according to the mathematical model of the first level.

Further improvement of the mathematical model of ischemic stroke will require the addition and solution of a system of differential equations, the number of which will correspond to the number of input variables. In such cases, you usually have to rely on approximate solutions.

To optimize the multistage scenario of preclinical trials of [²⁵Mg²⁺]₄PMC16, it is necessary to further develop the obtained model and combine it with the pharmacokinetic model of targeted delivery of [²⁵Mg²⁺]₄PMC16 NPs to the hypoxidated area of the brain at the border of healthy and necrotized tissues, for which it is necessary to refine and analyze the In Silico algorithms.

8. Combining mathematical models of individual processes into a system of differential equations

Let us detail the pathogenesis model of ischemic stroke and describe the hypothesis of the process at the cellular level in more detail. The previously used simplified model of ischemic stroke pathogenesis reduced death of brain cells from ischemia only to necrosis from hypoxia and did not take into account other factors of the process.

Let us modify the biological model as it was done in [11, 12] the model of brain cell (neurons, astrocytes and other glial cells) death, which takes into account both of the main mechanisms of cell death - apoptosis and necrosis. Nutrient deficiency promotes cell death programs initiation if the level of damage reaches a certain threshold value called D₀. At the present stage of research, we assume that with severe damage, apoptotic or necrotic variants of cell death can be realized with equal probability. The cells that die due to ischemia and subsequent intoxication initiate the immune system’s response. First of all, the body’s own defenders of nerve cells located in the brain, microglia, are activated, which can ingest (phagocytosis) the decay products, thereby participating in the removal of the destruction products and preventing the escalation of the inflammatory reaction. Activation of microglia prepares them for a state of readiness for phagocytosis and the synthesis of cytokines- specific proteins that coordinate the actions of cells of the immune system. Cytokines cause the accumulation of adhesion molecules in the vicinity of the damage, which promote the attachment (adhesion) of white blood cells traveling in the bloodstream to the endothelium of the blood vessel and their migration through the endothelium to brain cells. As a result, white blood cells overcome the hemato-encephalic barrier characteristic of a healthy body - a physiological mechanism that is designed to regulate the penetration into the brain of various substances introduced from the outside or circulating in the blood, in order to maintain the constancy of the physiological and physicochemical state of the brain. Thus, the cell adhesion molecules initiate the movement of white blood cells into the damaged area. It is known that neutrophils, the most numerous group of small white blood cells, as well as monocytes-macrophages, which are the largest white blood cells, play the main role in phagocytosis. Microglial cells, monocytes-macrophages and neutrophils contribute to the removal of dead apoptotic and necrotic brain cells from the body by phagocytosis. At the same time, activated microglial cells, pro-inflammatory cytokines and neutrophils themselves release toxic substances, which harm intact cells and contribute to the expansion of the brain lesion area [9].

Hence, the mathematical model of the dynamics of brain cell death can be supplemented with the corresponding equations and rewritten in the following form [11, 12]:

$$\begin{array}{c}\mathit{dH}/\mathit{dt}=-\mathrm{DH}\end{array}$$E18

$$\begin{array}{c}\mathrm{dN}/\mathrm{dt}=\left(1-p_A\right)\mathrm{DH}-\mathrm{\epsilon N}\end{array}$$E19

$$\begin{array}{c}\mathrm{d}{A}_s/\mathrm{dt}=p_A\mathrm{DH}-p_A\mathrm{D}\left(·-T_A\right)\mathrm{H}\left(·-T_A\right)\end{array}$$E20

$$\begin{array}{c}\mathrm{d}{A}_e/\mathrm{dt}=p_A\mathrm{D}\left(·-T_A\right)\mathrm{H}\left(·-T_A\right)-\mathrm{\epsilon }{A}_e\end{array}$$E21

$$\begin{array}{c}\mathrm{d}{M}_a/\mathrm{dt}=p{M}_a\left(c_A{A}_e+c_N\mathrm{N}\right)M_i-c_{\mathit{pro}}\left(M_a/T_{M,1}\right)\end{array}$$E22

$$\begin{array}{c}d{M}_i/\mathit{dt}=-p{M}_a\left(c_A{A}_e+c_{N}N\right)M_i+c_{\mathit{pro}}{M}_a{T}_{M,1}+\\ \left(c_{\mathit{Mi},1}{M}_i-c_{\mathit{Mi},2}{M}_{2i}\right)1_{t>\mathit{TM},2}-\epsilon M_i\end{array}$$E23

Where H — relative density of healthy brain cells, N — relative density of necrotic cells, A_s и A_e — relative density of brain cells that started and ended apoptosis, respectively, M_i и M_a — relative density of inactive and active microglial cells; the argument (· − T_A) represents a time delay T_A, corresponding to the characteristic duration of the cell apoptosis process; the sign 1_t>TM,2 also indicates the start of a number of important processes with a delay equal to T_M,2. In the equations of the mathematical model Eqs. (18)–(23), one of the key roles is played by the values of the specific rate of phagocytosis ε and the specific rate of cell death due to intoxication by decomposition products (inflammation)D, and the latter is different from zero only if the intoxication has reached the established limit level D₀:

$$\begin{array}{c}D=\left[(,p_{n,\left[\mathit{cy}\right]},C,+,p_{N,\mathit{Ln}},(,\left(L_n/\left(C_{L_n}+L_n\right)\right),\left(N+A_e\right),+,p_{N,N},N,),-,p_o,D_0\right]{}^{+}\end{array}$$E24

$$\begin{array}{c}\epsilon =e_{N,M_a}{M}_a+e_{N,L_m}{L}_m+e_{N,L_n}{L}_n+e_{N,M_i}{M}_i\end{array}$$E25

System (18–23) is supplemented by equations describing the dynamics of inflammation factors:

$$\begin{array}{c}d{L}_m/\mathit{dt}=c_{\mathit{Lm}}{M}_{\mathit{adh}}\left(·-T_{\mathit{Lm},\mathit{\text{in}}}\right)-p_{\mathit{dLm}}\left(L_m/T_{\mathit{Lm}}\right)\end{array}$$E26

$$\begin{array}{c}d{L}_n/\mathit{dt}=c_{\mathit{Ln}}{M}_{\mathit{adh}}\left(·-T_{\mathit{Ln},\mathit{\text{in}}}\right)-p_{\mathit{dLn}}\left(L_n/T_{\mathit{Ln}}\right)\end{array}$$E27

$$\begin{array}{c}\mathit{dC}/\mathit{dt}=\left(p_{\mathit{Ma},c}\left(M_a/M_a+c{M}_a\right)+{\mathit{pL}}_{m.c}\left(L_m/L_m+c{L}_m\right)\right)\left(N+A_e\right)-e_{\mathit{cy}}C\end{array}$$E28

$$\begin{array}{c}d{M}_{\mathit{adh}}/\mathit{dt}=\left[p_{\mathit{\text{Madh}},1}C-p_{\mathit{\text{Madh}},2}C{M}_{\mathit{adh}}-e_{\mathit{\text{Madh}}}{M}_{\mathit{adh}}\right]1_{\mathit{\text{vessel}}}\end{array}$$E29

Where L_m и L_n — the relative concentration of leukocytes of two types - monocytes-macrophages and neutrophils, respectively, C — relative concentration of cytokines, M_adℎ — the relative density of adhesion molecules. Eq. (29) describes the dynamics of cell adhesion molecules depending on the factors listed in its right-hand side only if the stroke nucleus is located in the vicinity of a blood vessel, and otherwise the density of adhesion molecules remains constant, corresponding to the initial condition. A feature of system (18–29) is the presence of functions with lagging arguments in the right-hand sides of the equations, where the lag T_Lm,n and T_Ln,in, as well as T_A, is due to biomedical considerations. The initial data for the components of the solution to the problem on the time interval t ∈ [−τ, 0] (where τ = max(T_A, T_Lm,in, T_Ln,in)), preceding the onset of the disease are set corresponding to the healthy state: the relative density of healthy cells H(t) = 1 and inactive microglia M_i(t) = 1, the values of the remaining variables of the problem are assumed to be zero; to simulate a stroke at t = 0 a given part of healthy cells (up to 40%) passes into necrotic and / or cells that have entered apoptosis.

Obviously, the mathematical description’s complexity will increase in proportion to the number of input variables and the increase in the number of differential equations. Simultaneously, the contribution of each of the new variables or a new differential equation describing the process under consideration at this stage is difficult to predict. In contrast, the entire set of processes accompanying ischemic stroke has not yet been sufficiently studied and described in the literature to be modeled in such detail.

9. Alternative ways of modeling ischemic stroke

At the same time, there are various models of population dynamics, which are used both in biology and in ecology and medicine while having sufficiently high reliability of the mathematical description of processes. Considering that any model is only a semblance of the original and the task of a complete repetition of a real object by a model is never set. Based on the existing models describing population dynamics, it seems possible to select the appropriate one and modify it for the tasks at hand. One of the variants of this approach is to use the modified Gompertz Equation [13], which describes the processes of population death rather well.

Hence, to optimize the multistage scenario of preclinical trials of [²⁵Mg²⁺]₄PMC16 it is necessary to combine the system of in silico stroke equations presented above with the pharmacokinetic model of targeted delivery of [²⁵Mg²⁺]₄PMC16 NPs to the ischemic penumbra zone.

10. Modeling the pharmacokinetics of fullerene-porphyrinic nano-cation exchangers carrying ²⁵Mg²⁺

The problems of modeling the pharmacokinetics of drugs have been studied extensively [14, 15, 16]. Such models are actively used in preclinical drug trials. Simultaneously, depending on the characteristics of the drugs under study, as well as the goals and objectives of such studies, one-chamber, two-chamber, three-chamber and four-chamber models are used.

The peculiarities of modeling the pharmacokinetics of fullerene-porphyrin nano-cation exchangers carrying ²⁵Mg²⁺ are that due to the spin effect of the ²⁵Mg²⁺ isotope, it hyperstimulates ATP synthesis, and due to the presence of the PMC16 “nanocontainer”, it has the property of “targeted” delivery to the area of the brain damaged by hypoxia.

In general, the dynamic processes of pharmacokinetics are modeled using systems of ordinary differential equations of the form [13]:

$$\begin{array}{c}\left\{\begin{array}{c}x=f\left(x\left(t\right),p\right)+{\sum }_{i=1}^{n}h\left(x\left(t\right),p\right)u_i\left(t\right)\\ y\left(t\right)=g\left(u\left(t\right),x\left(t\right),p\right)\end{array}\right.\end{array}$$E30

Where x(t) – n-dimensional function of the state (in pharmacokinetics - drug dose),

f(x(t), p) – a function that defines the structure of the model,

p – s-dimensional vector of parameters characterizing the process under consideration (in pharmacokinetics, the rate of drug transfer between organs),

h – a function that defines the structure of the input data,

u(t) – function of the input data (in pharmacokinetics - the method of introducing the drug into the body),

y(t) – k-dimensional function of experimental data (in pharmacokinetics - drug concentration in blood and/or urine),

g – a function that links the model to the dimensions.

By supplementing the system of equations Eqs. (18)–(29) with the system Eq. (30), we obtain a mathematical model of the In Silico I level of selective accumulation of cation-exchange PMC16 nanoparticles in brain cells and tissues for preclinical studies of the neuroprotective potential of fullerene-porphyrin nano-cation exchangers carrying ²⁵Mg²⁺ in relation to the pathogenesis of ischemic stroke.

11. Discussion of how this model will be correlated to real experiments

With regard to the correlation of the model with the data from real experiments, it is necessary to take into account a number of important circumstances arising from the specifics of the task. Namely: first of all, it is necessary to align the experimental results with the ischemic stroke model. This subproblem includes the coordination of each of the above-mentioned physiological processes that accompany the pathogenesis of this disease, expressed with separate differential equations. Furthermore, based on the in silico of stroke, which is consistent with the empirics, the same coordination of the equations describing the pharmacokinetics of PMC 16 is required.

However, this work must be completed. In this sense, further improvement of the model is planned in two main areas:

1. Complication of hypotheses used for the modeling processes and expansion of the system of differential equations;

2. Adaptation to problem conditions of existing semi-empirical models describing non-Markov population dynamics (Gomperz model, Verhulst logistic model, population size model in a periodic environment, population model with a smaller critical number, etc.) [17].

All these scenarios require coordination of in silico with experimental data obtained from in vivo of laboratory animals, which presupposes the following studies:

1. Defining the parameters of biological processes subject to experimental control in vivo.

2. Defining the variables and coefficients of the differential equations of the mathematical model, which are to be agreed with the experimental data.

3. Adaptation and optimization of relevant semi-empirical models (equations) of population dynamics and approximation of their parameters to the tasks set.

4. Determination of optimal mathematical methods for approximation and interpolation of experimental data.

5. Structure development and database creation for in silico goals and objectives.

6. Literature search and extraction of experimental data obtained by third-party researchers (external data). Their assortment, classification, and entry into specialized databases.

7. Comparison of external and internal (obtained as part of the framework of our own research) experimental data.

8. Output of transfer functions. Clarification of the pharmacokinetic equations in relation to the different methods of drug administration.

9. Preparation of algorithms for computer models.

This mathematical simulation (modeling) approach is an appropriate In Silico tool designed to describe and predict the key pharmacokinetic patterns of the in vivo distribution and the brain tissue accumulation of Magnesium-25 carrying — releasing PMC16 nanoparticles.

This tool seems promising for meeting the specific expectations of pharmacologists searching for the optimal, efficient and economical ways of planning this distinctive pharmacophore preclinical research.

12. Conclusions and prospects

This work presents algorithms for in silico modeling of selective accumulation of cation-exchange PMC16 nanoparticles in cells and tissues of the brain for preclinical studies of the neuroprotective potential of fullerene-porphyrin nano-cation exchangers carrying ²⁵Mg²⁺ concerning the pathogenesis of ischemic stroke. Solving this problem is extremely important for the optimization of multistage scenarios of preclinical trials of [²⁵Mg²⁺]₄PMC16 in experimental nanopharmacology of ischemic stroke. In the present study, these algorithms are in the spotlight.

As a result, we have obtained a relatively voluminous system of differential equations describing the pharmacokinetics of [²⁵Mg²⁺]₄PMC16 in relation to the pathogenesis of ischemic stroke.

In systems in which several processes are implemented simultaneously, the difficulty of accurately solving the modeling problem increases in proportion to the number of processes taken into account. The search for a solution to such systems by analytical methods is rather difficult. In practice, one usually has to rely on approximate solutions or the use of numerical methods and computer simulation technologies. For this reason, we did not search for an analytical solution to the presented system of differential equations in the framework of this work.

As an alternative way aimed at simplifying the system of differential equations underlying in silico, one can use combinations of slightly modified systems of Gompertz equations and the non-Markov concept of population dynamics.

In our previous works [1, 2, 3, 4, 5, 18], these algorithms were partially presented and are not described in detail in this study but are our research’s focus.

As seen from above, both Non-Markov population dynamics background and the Gompertz equation-based models were simultaneously applied here to harmonize a predicational pharmacokinetic validity and capabilities for the multivariant In Silico approach proposed for the Magnesium-25 releasing PMC16 nano-carriers as long as the latter are about to play a role of modulators of the brain hypoxia-related metabolic disorders.

Hereby, this is nothing more than an attempt to develop a simple and efficient computational tool applicable to optimize the decision-making process on exact steps and conditions of PMC16 engaging preclinical research strategy in the variable brain ischemia pharmacological studies.

Regarding the prospects for continuing work in this direction, it should be noted that at some point, as the algorithms become more complex, in silico formation without the use of artificial intelligence-based on computer neural networks will not be possible.

A further prospect of working on the In Silico project is to create conditions for bringing In Silico to the level of advanced smart technologies based on artificial intelligence neural networks.

Acknowledgments

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