Modular Model of Neuronal Activity That Captures the Dynamics of Main Molecular Targets of Antiepileptic Drugs

TL;DR

This paper introduces a modular model of neuronal activity that simulates how antiepileptic drugs affect key molecular targets in the brain.

Contribution

The novelty lies in integrating multiple molecular mechanisms into a single model to capture drug effects on epilepsy-related processes.

Findings

The model accurately simulates ion current and receptor dynamics relevant to epilepsy.

Simulations show consistent inhibitory effects on synaptic transmission, matching experimental data.

The model incorporates glutamate and GABA neurotransmitter interactions crucial for neural regulation.

Abstract

This paper presents a modular mathematical model of neuronal activity, designed to simulate the dynamics of main molecular targets of antiepileptic drugs and their pharmacological effects. The model was developed based on several existing synaptic transmission models that capture cellular processes crucial to the pathology of epilepsy. It incorporates the primary molecular mechanisms involved in regulating excitation and inhibition within the neural network. Special attention is given to the dynamics of ion currents (Na+, K+, Ca2+), receptors (AMPA, NMDA, GABAA, GABAB and mGlu), and neurotransmitters (glutamate and GABA). Examples of simulations illustrating the inhibitory effects on synaptic transmission are provided. The numerical results are consistent with experimental data reported in the literature.

Figures

Full text

1. Introduction

Epilepsy is a chronic neurological disorder characterized by a predisposition to spontaneous seizures, often accompanied by comorbidities such as anxiety, depression, and cognitive impairment [1]. In the brain of a person with epilepsy, excitatory activity dominates over inhibitory one, leading groups of neurons to synchronously produce powerful electrical discharges. These discharges can spread to other parts of the nervous system, triggering seizures. Seizures manifest as brief episodes of altered behavior, which can include symptoms such as loss of consciousness, convulsions, or sudden stiffening of the body [2].

The causes of epilepsy are diverse, including mutations, metabolic disorders, head injuries, infections, and other factors that impair brain function [3]. Treatment typically focuses on symptom suppression through antiepileptic drugs (AEDs; also termed antiseizure medications), but in severe cases, surgical intervention may be necessary to address the underlying causes of epileptic activity. However, the high toxicity of AEDs complicates the optimization of drug therapy strategies [4]. Despite the availability of over 30 antiepileptic drugs with various mechanisms of action, approximately one-third of patients are drug-resistant [5].

To enhance understanding of the mechanisms underlying epilepsy, mathematical models have been developed [6]. Such models are valuable tools for validating hypotheses, as they allow for detailed examination of neuronal dynamics under controlled conditions. Using computational modeling to study the mechanisms of neuronal response to epilepsy drug therapy can enhance general understanding of this complex disease, improve treatment efficacy, and reduce the risk of side effects.

There are three main approaches to brain modeling: nanoscopic (cellular), which focuses on individual neurons; microscopic (population), which models interactions among neuron populations; and mesoscopic (regional), which examines entire brain regions. The microscopic approach provides detailed insights into neural network dynamics [7,8]. However, simulating neuronal populations, which often involve tens of thousands of neurons, results in high computational complexity, limiting the practical applicability of this approach. Mesoscopic models can generate accurate electrophysiological brain signals, such as electroencephalogram signals, that are consistent with experimental recordings [9,10]. However, these models are often based on purely mathematical formulations, lacking direct biological interpretation, which complicates the validation of results. On the other hand, nanoscopic models simulate biophysical processes occurring within individual neurons, including receptor dynamics, intracellular signaling pathways, and molecular interaction networks, as well as their synaptic interactions, and therefore, in most cases, provide strong biological significance [11,12,13]. Although they are not capable of reproducing whole-brain electrophysiological signals, these models capture detailed dynamics of cellular processes and generate signals that can be experimentally validated in vitro.

In addition to classification by the level of system organization, mathematical models of neuronal activity can be broadly categorized into phenomenological and biophysical models. Phenomenological models describe neuronal dynamics using simplified mathematical representations with relatively few variables and parameters [9,10,14,15]. These models effectively capture general patterns of neural activity and can reproduce key features observed in experimental conditions. However, their abstract nature limits direct biological interpretation, which limits their applicability for studying the detailed mechanisms underlying molecular interactions. Biophysical models, in contrast, incorporate detailed descriptions of ion channel kinetics, neurotransmitter dynamics, and synaptic transmission [11,12,13,16,17,18,19]. By explicitly modeling the underlying cellular processes, these models provide a high degree of biological interpretability and allow for direct comparisons with experimental data. However, this increased biological realism comes at the cost of higher computational complexity and a greater number of parameters. Typically, phenomenological models are used in mesoscopic and microscopic approaches, whereas biophysical models are predominantly applied at the nanoscopic level.

Consequently, nanoscopic biophysical models are particularly well-suited for simulating the dynamics of molecular targets of AEDs, which is the primary focus of this study.

Computational biology methods have developed very rapidly in recent years, and breakthroughs in neuronal activity modeling have enabled direct correspondences between experimentally measurable quantities and simulated biophysical processes, such as changes in ion concentrations and membrane potential oscillations [20,21,22]. This capability facilitates biologically meaningful modeling of the biophysical pathways targeted by AEDs, thereby simplifying the optimization of drug therapy.

The objective of this study was to develop a mathematical model of neuronal activity that accurately describes the functioning of neurons, taking into account the dynamics of the main molecular targets of AEDs. To accomplish this, we reproduced results from several existing synaptic transmission models on the BioUML (Biological Universal Modeling Language) platform, focusing on key cellular processes crucial in the pathology of epilepsy. Based on reproduced models, we constructed a modular model.

The paper is organized as follows. Section 2 presents the main results of the study. It first characterizes the main molecular targets of antiepileptic drugs, whose dynamics were considered in the model’s development, and formulates the corresponding requirements for the model. This is followed by a detailed description of the modular structure of the model and its main equations. The section then presents the results of numerical simulations, including the dynamics of neuronal and astrocytic activity, ion and neurotransmitter concentrations, demonstration of interneuron-mediated inhibition, and an example of drug effect simulation, with comparisons to experimental data from the literature. Section 3 discusses the obtained results and outlines the limitations and potential extensions of the model. Section 4 briefly describes the features of the BioUML software used for building the modular model and performing numerical calculations. Section 5 summarizes the main conclusions of the study. Appendix A, Appendix B, Appendix C, Appendix D and Appendix E provide a detailed description of all the model equations, as well as the parameter values and initial conditions for the variables.

2. Results

2.1. Main Molecular Targets of Antiepileptic Drugs

The first-line treatment for epilepsy typically involves AEDs, which symptomatically suppress seizures. Their mechanisms of action are aimed at regulating synaptic activity, modifying the intrinsic excitability properties of neurons and balancing inhibitory and excitatory neurotransmission, thereby reducing the probability of seizure occurrence.

The actions of antiepileptic drugs are based on the effect on molecular targets that regulate neuronal activity. Main molecular targets of AEDs include voltage-gated sodium (Na^+^), potassium (K^+^) and calcium (Ca^2+^) channels; ionotropic receptors of a-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA-Rs), N-methyl-D-aspartate (NMDA-Rs) and γ-aminobutyric acid (GABA_A_-Rs); as well as synaptic vesicle protein 2A (SV2A), GABA transporter 1 (GAT-1) and GABA aminotransferase (GABA-T) [23]. These targets are schematically shown in Figure 1.

In addition to GAT-1, which is predominantly expressed in neurons and acts as the main regulator of extracellular GABA in the hippocampus, there is also GAT-3, primarily expressed in glial cells, which plays a significant role in conditions of excessive neuronal activity [24]. The experimental observations indicate GAT-3 can release GABA from astrocytes into the synaptic space, contributing to astrocytic modulation of hyperexcitability [25,26]. Although no current antiepileptic drugs specifically target GAT-3, its activation leads to the inhibition of synaptic transmission [24]. Thus, GAT-3 dynamics are also of significant interest in modeling excessive synaptic activity.

Pharmacological action on primary targets of AEDs suppresses seizures through the following mechanisms:

• Activation of voltage-gated K^+^ channels or inhibition of voltage-gated Na^+^ and Ca^2+^ channels, leading to membrane hyperpolarization and thereby returning the neuron to its resting state;
• Enhancement of GABA-mediated inhibition via inhibition of the GAT-1 GABA transporter in neurons and astrocytes, inhibition of GABA transaminase, and activation of GABA_A_ receptors;
• Attenuation of synaptic excitation through inhibition of postsynaptic AMPA-Rs and NMDA-Rs;
• Direct modulation of synaptic release by activating presynaptic SV2A and inhibiting voltage-gated Ca^2+^ channels.

2.2. Requirements for the Modular Model

Based on the mechanisms of action of AEDs discussed above, the developed modular model of neuronal activity should account for the following:

• Na^+^, K^+^ and Ca^2+^ concentration dynamics, which are fundamental to the generation of action potentials (APs). Specifically, the model should incorporate voltage-gated channels for these ions.
• Dynamics of AMPA-Rs, NMDA-Rs and GABA-Rs, which are key regulators of excitatory and inhibitory synaptic transmission.
• Astrocytic regulation of synaptic cleft homeostasis, particularly the dynamics of the GAT-3.

These aspects have been incorporated into the model’s design.

2.3. Modular Model Structure

It is well established that glial cells, particularly astrocytes, play a crucial role in neurotransmission in addition to neurons [27]. Therefore, to accurately describe the mechanisms of neuronal activity, it is appropriate to consider tripartite (three-part) synapse models, which include not only the presynaptic and postsynaptic neurons but also perisynaptic astrocytes.

In our analysis of existing neuronal activity models, we focused on those within the tripartite synapse formalism.

Table 1 lists the neuronal activity models that were most relevant for our study and directly used in the development of the modular model. While numerous other models have been proposed in the literature, offering valuable insights into various aspects of neuronal dynamics, they were less suitable for direct integration into our framework. A broader discussion of these models and their relevance to our specific objectives is provided in Section 3.

The mathematical model of neuronal activity that captures the dynamics of main molecular targets of antiepileptic drugs was developed based on the principles of modular modeling of complex biological systems using the BioUML platform. The modular structure of the model in BioUML is shown in Figure 2.

The core of the model was based on [11] and was further extended to include the regulation of GABA and glutamate dynamics by astrocytes as detailed in [25]. Additionally, the model incorporated GABA_A_, GABA_B_, and metabotropic glutamate receptors (mGlu-Rs) [25,28], as well as the inhibitory effects of interneurons [12].

The BioUML implementation of the model consists of six modules: presynaptic and postsynaptic neurons, perisynaptic astrocyte, interneuron, extracellular space, and externally applied current. It includes 247 variables, 41 differential equations, 174 parameters, and 235 assignment operations. Hereafter, it is referred to as the modular model of neuronal activity.

The model’s modules are briefly described as follows:

1. Presynaptic neuron (Presynaptic_neuron): Module that describes the dynamics of the presynaptic axon terminal, incorporating voltage-gated ion channels, sodium-potassium and calcium ATPases (NKA and PMCA), GABA_A_-Rs involved in tonic astrocyte-mediated modulation, and external stimulus. It calculates the concentration of glutamate released into the extracellular space, with the release rate primarily determined by presynaptic calcium entry, while being modulated by the dynamics of mGlu-Rs and GABA_B_-Rs. Ion fluxes across the membrane are calculated as well.
2. Postsynaptic neuron (Postsynaptic_neuron): This module describes the dynamics of the postsynaptic dendritic spine, including voltage-gated ion channels, NKA, PMCA, AMPA-Rs, NMDA-Rs, as well as GABA_A_-Rs activated by astrocyte-mediated tonic GABA signaling, and a weak external stimulus. Ion fluxes across the membrane are also computed.
3. Perisynaptic astrocyte (Astrocyte): This module describes the dynamics of the astrocyte, including the excitatory amino acid transporter 2 (EAAT-2), the sodium-calcium exchanger (NCX), NKA, inwardly rectifying potassium channel (Kir4.1), and the GABA transporter 3 (GAT-3). Ion diffusion within the astrocytic process is considered, along with ion and neurotransmitter fluxes across the membrane.
4. Interneuron (Presynaptic_interneuron): This module calculates the inhibitory effect of the interneuron by calculating the release of GABA into the extracellular space. The excitation is initiated by an external current injection and mediated by the activation of voltage-gated ion channels.
5. Extracellular space (Extracellular_space): This module describes the dynamics of the synaptic cleft, where concentrations of ions and neurotransmitters are established based on the fluxes across the membranes of the adjacent cellular compartments. Ion diffusion is also considered.
6. Externally applied electrical current (Applied_current): This auxiliary module calculates stimulating currents for the presynaptic neuron, postsynaptic neuron, and interneuron.

2.4. Main Equations of the Modular Model

This section provides a brief overview of the main equations underlying the modular model of neuronal activity (Equations (1)–(24)). Detailed descriptions of all equations, parameter values, and initial conditions for the variables can be found in Appendix A, Appendix B, Appendix C, Appendix D and Appendix E.

The presynaptic and postsynaptic neurons, along with the interneuron, are modeled using the Hodgkin-Huxley-based description [29] for voltage-gated ion channels dynamics. The membrane potential of these cells is calculated in Equations (1)–(3) as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ , $[eqn]$ and $[eqn]$ denote the membrane potentials of the presynaptic neuron, postsynaptic neuron and interneuron, respectively; $[eqn]$ is the membrane capacitance; $[eqn]$ , $[eqn]$ and $[eqn]$ are the currents of Na^+^, K^+^ and Ca^2+^ ions across the membrane, respectively; $[eqn]$ is the leak current; $[eqn]$ and $[eqn]$ are the externally applied currents. $[eqn]$ is used to excite the presynaptic neuron and interneuron, while $[eqn]$ is required to remove the NMDA-Rs magnesium block. $[eqn]$ is introduced following the original synapse-astrocyte model [11] as a numerical surrogate for local postsynaptic depolarization in the absence of explicit somatic and dendritic morphology. This weak depolarizing current is applied with a 2 ms delay relative to the presynaptic neuron and has a much smaller amplitude, calculated to ensure that the postsynaptic neuron is excited as a result of presynaptic excitation, but that $[eqn]$ alone is insufficient to trigger postsynaptic neuron’s excitation.

Presynaptic membrane ionic currents are calculated considering voltage-gated ion channels, NKA, PMCA, and leak channels using Equations (4)–(6):

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ , $[eqn]$ , $[eqn]$ denote voltage-gated Na^+^, K^+^ and Ca^2+^ channels currents, respectively; $[eqn]$ , $[eqn]$ are NKA and PMCA generated currents; and $[eqn]$ , $[eqn]$ , $[eqn]$ represent leak currents.

The presynaptic neuron also calculates the concentration of glutamate released into the synaptic cleft (Equation (7)), with its release probability primarily determined by presynaptic calcium entry and modulated by activation of mGlu-Rs and GABA_B_-Rs (Equation (8)), where the former increases [30] and the latter decreases it [31].

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ denotes the concentration of glutamate released from the presynaptic neuron vesicles per second; $[eqn]$ and $[eqn]$ are the probabilities of glutamate release at rest and over time, respectively; $[eqn]$ is glutamate release rate; $[eqn]$ is glutamate concentration per vesicle; $[eqn]$ and $[eqn]$ are the fractions of activated mGlu-Rs and GABA_B_-Rs, respectively.

Ionic currents of the postsynaptic neuron (Equations (9)–(11)) additionally include AMPA-Rs and NMDA-Rs, which are necessary for successful synaptic transmission.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ and $[eqn]$ are AMPA-Rs and NMDA-Rs generated currents, respectively.

The interneuron inhibits synaptic transmission by releasing GABA into the extracellular space, with its released concentration described by Equation (12). This process is modeled using the Tsodyks–Markram synaptic plasticity model [32], which accounts for the limited availability of synaptic resources and their recovery dynamics (Equations (13)–(15)).

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Within this scheme $[eqn]$ denotes the concentration of GABA released from the interneuron per second; $[eqn]$ , $[eqn]$ and $[eqn]$ are the fractions of synaptic vesicles in the recovered, active and inactive states, respectively; $[eqn]$ and $[eqn]$ are vesicle recovery and inactivation time constants, respectively; $[eqn]$ is the number of docked vesicles; $[eqn]$ is the GABA concentration in single vesicle; $[eqn]$ is the Dirac delta function; $[eqn]$ —interneuron’s spike times.

Astrocytes play a crucial role in maintaining the homeostasis of ions and neurotransmitters at the synapse [33]. Each astrocyte extends several primary processes (defined as branches), which split into increasingly smaller processes, eventually ending in peripheral processes with sheet-like structures known as leaflets. These leaflets lack organelles but have various transporters dwelling in their membrane [34,35]. Astrocytic leaflets connect with synapses thus forming a tripartite synapse and maintaining synaptic microenvironment by regulating ionic and neurotransmitters concentrations.

In this study, the perisynaptic cradle (PsC) refers to the astrocytic compartment that surrounds synapses and regulates ion and neurotransmitter concentrations in the synaptic cleft. The PsC includes an astrocytic leaflet as well.

The ion and neurotransmitter currents within the PsC are calculated using Equations (16)–(20):

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ denote EAAT-2, NCX, Kir4.1 and GAT-3 generated currents, respectively.

Ion diffusion within the leaflet is also considered, as described by Equation (21):

$$[eqn]$$

Here, $[eqn]$ denotes leaflet diffusion current; $[eqn]$ represents equilibrium potentials for ions relative to the intracellular space of the astrocyte; $[eqn]$ is the leaflet length; $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ denote the Poole–Frankel channel constant, well activation energy, dynamic permittivity, elementary charge, Boltzmann’s constant and absolute temperature, respectively.

The dynamics of the extracellular space (ECS) are determined by considering ion and neurotransmitter currents from all other compartments, as well as ion diffusion within the cleft, as described by Equations (22) and (23):

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ represents ECS diffusion current and $[eqn]$ denotes the concentration of neurotransmitters released by the presynaptic neuron and interneuron.

ECS diffusion is calculated relative to the global extracellular space (GECS), which refers to the extensive extracellular region surrounding neurons and glial cells, as described by Equation (24):

$$[eqn]$$

where $[eqn]$ denotes equilibrium potentials for ions relative to the GESC; $[eqn]$ is the leaflet cross-sectional area; $[eqn]$ and $[eqn]$ denote conductance scaling factor and maximal diffusion conductance, respectively. It is assumed that the GESC has a much larger volume than the ECS; therefore, the concentrations of all ions in the GECS are assumed to remain effectively constant.

2.5. Numerical Simulation Results

2.5.1. Preliminary Information

The results of the numerical simulations for the modular model of neuronal activity are presented for the entire simulation period of 100 s. The presynaptic and postsynaptic neurons are stimulated by an externally applied pulsed current for 50 s, starting at the 10 s. Inset graphs display the last second of stimulation to provide a more detailed view of the dynamics of currents, concentrations, and membrane potentials. Positive values of transmembrane currents indicate exiting the cellular compartment (outward flux), while negative values represent entering the cellular compartment (inward flux).

2.5.2. Neuronal Activity

The presynaptic neuron (Pre) is stimulated with a pulsed depolarizing current of 10 µA/cm^2^, at a frequency of 10 Hz, and a pulse width of 3 ms. This frequency was chosen to provide stable synaptic activation and is commonly used in computational modeling studies of tripartite synapse dynamics [11,25,36]. The stimulation causes membrane depolarization (Figure 3A), which triggers the opening of voltage-gated Na^+^, K^+^ and Ca^2+^ channels, leading the neuron to generate action potentials. The action potential of the Pre has a magnitude of approximately 110 mV, with ~100 mV depolarization and ~10 mV afterhyperpolarization.

Glutamate and GABA are the brain’s primary excitatory and inhibitory neurotransmitters, respectively [37]. During each action potential, Ca^2+^ enters the Pre through voltage-gated Ca^2+^ channels (Figure 3B), leading to glutamate release into the synaptic cleft. Glutamate then activates AMPA-Rs and NMDA-Rs on the postsynaptic neuron (Post). Activation of AMPA-Rs allows Na^+^ to flow into the Post (Figure 3C), while NMDA-Rs permit both Na^+^ and Ca^2+^ to enter the cell (Figure 3D), with K^+^ efflux into the synaptic cleft (Figure 3E). The activity of glutamate receptors depolarizes the Post membrane, triggers the opening of voltage-gated ion channels, and initiates action potentials (Figure 3F). In contrast, GABA serves as an inhibitory neurotransmitter by binding to GABA_A_-Rs and GABA_B_-Rs, thereby retarding synaptic transmission. Under normal conditions, GABA levels in the ECS remain low, providing a mild but consistent inhibitory effect on nearby neurons. However, external sources of GABA, such as interneurons, can significantly increase its concentration, leading to a substantial decrease in synaptic activity (as shown in Section 2.5.5).

To remove the magnesium block of NMDA-Rs, the Post is stimulated with a weak depolarizing current simultaneously with glutamate release from the Pre. This depolarization alone is insufficient to induce postsynaptic firing and does not substitute synaptic excitation mediated by AMPA-Rs activation. Notably, the slight reduction in action potential amplitudes in both the Pre and Post after peaking at the onset of stimulation may reflect the mechanism of short-term depression [38].

2.5.3. Astrocytic Regulation

Figure 4 shows the ion and neurotransmitter currents through the main astrocytic pathways: EAAT-2, NCX, NKA, Kir4.1 and GAT-3. Glutamate released into the ECS during Pre stimulation is taken by the astrocytic EAAT-2 (Figure 4A), leading to an influx of Na^+^ (Figure 4B) and an efflux of K^+^ (Figure 4C). According to published data [25], the Na^+^ influx through EAAT-2 is sufficient to activate the GAT-3, which releases GABA (Figure 4D) and Na^+^ (Figure 4E) into the synapse, contributing to neuronal inhibition. The maintenance of ionic homeostasis also involves NKA (Figure 4F,G), NCX (Figure 4H,I) and Kir4.1 (Figure 4J). The current through NCX is significantly lower than that through EAAT-2, consistent with experimental observations [39]. The slow release of GABA by the GAT-3 contributes to long-lasting tonic inhibition of nearby neurons, sustaining a low yet steady GABA level in the synaptic cleft [26,40,41]. This contrasts with the transient, or phasic, inhibition typically caused by exocytosis [41].

2.5.4. Ion and Neurotransmitter Concentrations

Figure 5 illustrates the dynamics of ionic concentrations in the presynaptic and postsynaptic neurons. Resting ion and neurotransmitter concentrations were set according to experimentally derived physiological data reported in the literature [42,43]. Membrane depolarization causes the opening of voltage-gated ion channels, allowing ion currents to flow according to their concentration gradients, resulting in Na^+^ and Ca^2+^ influxes and K^+^ efflux (Figure 5A–F). Once the stimulation ends, the ion concentrations quickly return to their equilibrium levels due to the activity of NKA and PMCA.

The dynamics of ion and neurotransmitter concentrations in the PsC and ECS are shown in Figure 6. Synaptic activity alters synaptic concentrations, with astrocytes actively regulating synaptic homeostasis.

When presynaptic vesicles release their contents, glutamate concentrations in the synaptic cleft quickly peak around 1 mM (Figure 6H). This glutamate is rapidly cleared from the synapse by the astrocytic EAAT-2 and subsequently decays (Figure 6G). The clearance time for glutamate from the synaptic space is approximately 30 ms (Figure 6H, lower inset), consistent with experimental observations [44]. During each EAAT-2 transport cycle, Na^+^ is transported into the PsC (Figure 6A), while K^+^ is moved out to the ECS (Figure 6D). Astrocytic NKA prevents excessive K^+^ accumulation in the ECS, which could result from NMDA-Rs, EAAT-2 and Kir4.1 (Figure 6C), and facilitates Na^+^ efflux from the astrocyte (Figure 6B). Changes in extracellular K^+^ concentration remain within the acceptable increase range of several millimolars during periods of intense neuronal activity [45]. The activity of NCX is highly dependent on the intra- and extracellular Na^+^ concentrations, so the dynamics of Ca^2+^ concentrations are mainly determined by Na^+^ levels in the PsC and ECS (Figure 6E,F).

In the absence of inhibitory input from interneurons, extracellular GABA concentration is maintained at a low but constant level (around 1 µM, Figure 6J). The weak current of GABA is not sufficient to significantly change its concentration in the PsC (Figure 6I). The results of the simulation under conditions of interneuron inhibition are presented in the following section.

After stimulation ends, all concentrations return to their equilibrium levels. Na^+^ and Ca^2+^ concentrations quickly recover due to the activity of NKA and NCX, which pump these ions out of the PsC. K^+^ concentration recovery is slower, as K^+^ continues to shuttle between the ECS and PsC through NKA, Kir4.1 and leak channels.

2.5.5. Interneuron Inhibition

At rest, GABA concentration in the synaptic cleft remains low (around 0.8 µM [46]). Interneurons are able to provide inhibition to neuronal networks and dictate the temporal pattern of activity of pyramidal neurons [46]. Stimulation of the interneuron leads to the opening of voltage-gated channels and the subsequent exocytosis of synaptic vesicles containing GABA. This results in a rapid increase in GABA concentration in the synaptic cleft, reaching levels far above the baseline (tens of µM, Figure 7A). The large fluctuations in GABA concentration are driven by the active uptake mechanism of the GAT-3 transporter.

In the ECS, GABA binds to GABA_A_-Rs and GABA_B_-Rs on neuronal surfaces, generating inhibitory currents and reducing the likelihood of glutamate release. These factors lead to a decrease in the frequency and amplitude of membrane potential oscillations (Figure 7B,D) and a reduction in glutamate synaptic concentration (Figure 7C). Gradually, the membranes of both neurons stop depolarizing beyond the threshold, halting action potential generation and thereby stopping synaptic transmission.

The Pre remained in an excited state longer than the Post because it is affected by an external stimulus with high amplitude, whereas the only source of excitation of the Post is the activation of glutamate receptors.

It should be noted that the accumulation of excessive GABA concentration in the synapse can also be caused by a number of other factors, including the inhibition of GABA transporter GAT-1 and GABA aminotransferase GABA-T, which are the basis for the mechanisms of action of some antiepileptic drugs (e.g., valproic acid [47]). While these components are not considered in our model, we assume that the increase in GABA levels through interneuron release has a similar impact on system dynamics, as the primary outcome of targeting GAT-1 and GABA-T is the elevation of GABA concentration in the synaptic cleft.

2.5.6. Simulation of Drug Effect

One of the most common targets of AEDs is voltage-gated Na^+^ channels [48]. Inhibitors of these channels include drugs like carbamazepine, lamotrigine, phenytoin, lacosamide, and others [23]. These drugs bind predominantly to depolarized Na^+^ channels, shifting them into a non-conducting state similar to inactivation but with slower recovery. Inhibitory effect occurs during prolonged or repetitive neuronal depolarization, which is crucial during epileptic seizures [48]. The drug concentration is calibrated to prevent pathological activity without completely blocking neuronal function. This approach ensures a balance between therapeutic effectiveness and the preservation of normal neuronal activity.

To model the effect of Na^+^ channel inhibitors, one can reduce the conductance of these channels and/or adjust the parameters governing their activation and inactivation kinetics. As an example, we will gradually decrease the maximal conductance of voltage-gated Na^+^ channels starting from the point of stimulus application, thereby simulating the selective enhancement of their slow inactivation.

Figure 8 shows the simulation results under conditions of voltage-gated Na^+^ channels inhibition. When the maximal conductance is reduced by less than 20% (Figure 8A), the Na^+^ channels still provide enough current to generate action potentials (Figure 8B). With a reduction of more than 20% (Figure 8A), the Pre continues to generate action potentials in the presence of a strong external stimulus (Figure 8B), but the glutamate concentration released (Figure 8C) is insufficient to excite the Post (Figure 8D). Reduction of more than 30% (Figure 8A) leads to complete inhibition of synaptic transmission (Figure 8D).

A study on differentiated NG108-15 cells, which are frequently used in in vitro studies instead of primary-cultured neurons, demonstrated that a decrease in voltage-gated Na^+^ channel availability suppresses action potential generation, while a sufficient level of Na^+^ conductance is required to sustain neuronal excitability [49]. In our model, a gradual reduction in maximal Na^+^ channel conductance leads to a progressive decrease in action potential amplitude. Once a critical threshold (~30% reduction) is reached, action potential generation ceases entirely, as the membrane potential no longer surpasses the firing threshold. This result supports the idea that Na^+^ channel inhibition can significantly impact neuronal firing and synaptic function.

Thus, the model is capable of simulating the basic effects of pharmacological treatments. Given the complexity of AED mechanisms and the limited availability of experimental data, a thorough investigation requires extensive parameter optimization and validation. This aspect will be the focus of our future work. Meanwhile, the developed modular model serves as a valuable standalone framework for simulating neuronal activity while accounting for the dynamics of the main molecular targets of AEDs.

3. Discussion

Synaptic transmission is a widely researched area in both experimental and computational research. In this study, we developed a mathematical model of neuronal activity that incorporates the dynamics of main molecular targets of AEDs, the role of astrocytic transporters in regulating ion and neurotransmitter concentrations, and the impact of interneuron inhibition.

The constructed model describes synaptic activity by calculating changes in ion and neurotransmitter concentrations which are consistent with experimental data. A key advantage of the model is that it includes the dynamics of main molecular targets of AEDs, which opens up opportunities for modeling drug effects on synaptic transmission, ranging from changes in ion and neurotransmitter concentrations to broader impacts on neuronal networks. Specifically, the model explicitly accounts for voltage-gated Na^+^, K^+^ and Ca^2+^ channels, AMPA, NMDA, GABA_A_ and GABA_B_ receptors, as well as key astrocytic transporters, including the glutamate transporter EAAT-2 and the GABA transporter GAT-3.

For model construction, we selected biophysically detailed nanoscopic-level neuronal activity models that incorporate the dynamics of the main molecular targets of AEDs and were best suited for integration into a unified framework. These selected models are listed in Table 1.

Beyond the models integrated into the modular model, we acknowledge the existence of many other noteworthy neuronal activity models. However, they were not directly suitable for integration due to specific limitations. For example, some models focus exclusively on the detailed dynamics of a single neuron [13,50,51], or describe astrocytic regulation of presynaptic activity [52,53,54] without considering synaptic transmission, which is crucial in the context of epilepsy. Additionally, some models adopt a population-based approach, averaging neuronal ensemble activity instead of capturing the dynamics of individual neurons [55,56]. While these models provide valuable insights into specific aspects of neuronal dynamics, our approach aims for a detailed representation of synaptic transmission at the level of single neurons. Finally, several biophysical models of the tripartite synapse have been proposed [11,57,58,59], which could have been used as a foundation for our modular model. However, none of them simultaneously account for all the main molecular targets of AEDs considered in our study. We selected [11] as the most suitable basis because it already incorporated many of the key processes relevant to our objectives while omitting unrelated mechanisms. This choice minimized the need for additional modifications to incorporate the main molecular targets of AEDs while preserving biological accuracy.

The current formulation of the model incorporates a relatively simple Hodgkin-Huxley formalism to describe the dynamics of voltage-gated ion channels. This approach provides a biologically meaningful representation of action potential generation and ionic currents, which is sufficient for capturing the dynamics of the main molecular targets of antiepileptic drugs. However, this approach has inherent limitations in reproducing finer effects of pharmacological interventions, such as representation of drug-binding and unbinding kinetics, voltage-dependent channel affinity, or plasma concentration-time profiles. The goal of the present study is not to reproduce detailed experimental or clinical pharmacokinetic/pharmacodynamic relationships, but to establish a biophysically detailed model that explicitly accounts for the dynamics of all the main molecular targets of AEDs, including ion channels, synaptic receptors, and astrocytic transporters. This foundation is necessary to enable systematic extension and refinement of the model in future studies.

Similarly, the model employs simplified representations of GABA_B_ and mGlu receptor activation, which limits its ability to capture the full diversity of downstream signaling pathways associated with these receptors. Clarification and explicit modeling of the corresponding molecular cascades constitute an additional direction for further development of the model.

Also, at the present stage, the model does not explicitly reproduce epileptiform discharges or ictal–interictal transitions; its primary purpose is to capture the dynamics of the main molecular targets of antiepileptic drugs under physiologically relevant neuronal activity. In the future, we plan to add a pathological excitation pathway to the model with the aim of simulating seizures and exploring epilepsy drug therapy. The extended model will then be integrated into our previously published multilevel model of epileptic seizures [60] to study drug therapy effects across different levels of brain organization.

4. Methods

4.1. Modeling Software

The construction and numerical calculations of the model were carried out on the BioUML platform (https://sirius-web.org/bioumlweb/, version 2025.2, accessed on 1 January 2026), an open-source Java-based integrated environment designed for modeling complex biological systems and analyzing biomedical data [61,62]. In this section, we provide a brief overview of the specific functionality of the BioUML platform that was used in the creation of the model.

4.2. Visual Modeling

Graphical representation simplifies the development of mathematical models for biological systems. In the visual modeling approach, processes and systems are depicted as diagrams, which serve as the basis for numerical simulations.

The model was constructed on the BioUML platform using visual modeling as follows. Initially, the biological system is represented as a diagram, from which BioUML generates Java code that describes the corresponding system of differential equations. This code is then compiled and executed by the numerical solvers integrated into BioUML to simulate the system’s dynamics.

4.3. Modular Modeling

The model was developed using a modular approach implemented in the BioUML platform, which allows any complex system to be broken down into elementary components called modules (also referred to as submodels or blocks). Each module can be created and configured independently, using its own formalism and level of detail. The interaction between modules is defined by the connections between their internal variables. These connections represent signal transmission pathways between the corresponding blocks. The graphical elements of the modular approach used in model construction are presented in Table 2.

Constructing a modular model involves defining submodels, specifying their input and output variables using ports, and connecting these ports through appropriate connections. Each module can have multiple connections with other modules. Additionally, in order to simplify the visual representation of the model, a bus element was used in the modular diagram to allow remote connections between modules.

4.4. Numerical Solution

One unit of model time corresponds to one second of real time. All simulation results presented in this work were performed using a forward Euler numerical integration scheme with a fixed time step $[eqn]$ (1 μs), with the simulation time set to 100 s.

5. Conclusions

Within this paper, the modular computer model of neuronal activity has been presented. It captures the dynamics of the main molecular targets of antiepileptic drugs, including voltage-gated Na^+^, K^+^ and Ca^2+^ channels, as well as AMPA, NMDA and GABA receptors.

The model simulates the signal transmission from the presynaptic to the postsynaptic neuron across the synaptic cleft, where the dynamics of ion and neurotransmitter concentrations are actively regulated by astrocytes through various channels and transporters, including the GABA transporter GAT-3, which is particularly important in modeling excessive synaptic activity. Additionally, the system incorporates an interneuron that can inhibit synaptic transmission by releasing GABA.

One of a key advantage of the model is its ability to simulate the effects of antiepileptic drugs on their molecular targets. As an example, the results of simulation under conditions of voltage-gated Na^+^ channels inhibition were presented.

Comparative analysis of the model’s numerical results with data from the literature revealed realistic reproduction of ion and neurotransmitter concentration dynamics, as well as currents through astrocytic channels, transporters and exchangers. The developed model is well-suited for use in the study of epilepsy drug therapy.