We represent the electrical potential of neurons and their surrounding extracellular space using the Extracellular-Membrane-Intracellular (EMI) model (Agudelo-Toro and Neef, 2013; Tveito et al., 2017b; Jæger and Tveito, 2020).

The EMI model is given by the following system of equations:

Here, u_i and u_e (in mV) are the intracellular and extracellular potentials, defined in the intracellular and extracellular domains, Ω_i and Ω_e, respectively. Moreover, v = u_i−u_e is the membrane potential given on the membrane, Γ, which is defined as the interface between Ω_i and Ω_e. In addition, σ_i and σ_e (in mS/cm) are the intracellular and extracellular conductivities, respectively, n_i and n_e are the outward pointing normal vectors of Ω_i and Ω_e, respectively, C_m (in μF/cm²) is the specific membrane capacitance, and I_s (in μ A/cm²) is a stimulation current density. Unless otherwise stated, I_s = 0. Time, t, is given in milliseconds (ms). Furthermore, I_ion (in μ A/cm²) is the current density through different types of ion channels on the cell membrane, s are gating variables and ionic concentrations involved in modeling these current densities, and F(t, v, s) model the dynamics of these additional state variables.

For the model neurons considered in this study, the density of the different types of ion channels vary for different locations of the neuron, like the soma, the dendrites and the axon initial segment (AIS). This results in a spatially varying expression for I_ion, which is described in more detail in the following subsections and in Supplementary material.

In addition, we consider axons partially covered by myelin. In the myelin covered part of the membrane, both I_ion and C_m are set to zero, resulting in the no-flux boundary conditions

on the myelinated membrane.

In our simulations, we let the neurons be surrounded by approximately 500 μm of extracellular space in spatial each direction. On the outer boundary of the extracellular domain, we apply the Dirichlet boundary condition

We model cerebellar Purkinje neurons using the EMI modeling setup described in Jæger and Tveito (2025b), based on the cable model formulation from Masoli et al. (2024). The Purkinje neuron is divided into seven different morphological parts, the dendrites, the soma, the axon initial segment (AIS), a ParaAIS, three Ranvier nodes separated by myelinated membrane sections and an axon collateral (see Figure 1A). The ion channel densities for each of these sections are described in Supplementary material. The three-dimensional (3D) finite element mesh of a neuron and the surrounding extracellular space is illustrated in Figure 1B.

In the absence of stimulation, the Purkinje neuron fires spontaneous APs. These are initiated in the AIS and propagate left through the soma and dendrites and right through the axon (see Figure 1E). In the soma and AIS clear negative EPs are generated during AP firing, whereas positive EPs are generated near the dendrites. Figure 1D provide spatial solution snapshots of the membrane potential and the surrounding EP at three points in time during AP firing.

If a positive stimulation current, I_s, is applied in the soma of the Purkinje neuron model, the frequency of the AP firing is reduced, and if the stimulation is sufficiently strong (e.g., I_s = 0.5 μ A/cm²), the AP firing is completely suppressed (see Figure 1C).

The model for neocortical layer 5 pyramidal neurons is set up in a similar manner as the model for cerebellar Purkinje neurons (see Figure 2). The model for the membrane dynamics is based on Hay et al. (2011), and the specific parameter values are specified in Supplementary material.

The pyramidal neuron morphology consists of eight different parts, the apical dendrites, the basal dendrites, the soma, the axon neck between the soma and the AIS, the AIS, three Ranvier nodes separated by myelinated membrane sections and an axon collateral (see Figure 2A). The channel densities and geometry for each of these different sections are specified in Supplementary material. In Figure 2B, the 3D finite element mesh used to represent a pyramidal neuron and its surrounding extracellular space is illustrated.

In the absence of stimulation, the pyramidal neuron do not fire APs, but when the soma is stimulated by a sufficiently strong negative stimulation current, I_s, APs are fired (see Figure 2C). When the strength of I_s is further increased, the frequency of the AP firing increases.

In response to a somatic stimulation (I_s = −27 μA/cm²), the AP is first initiated in the AIS and propagates left through the axon neck, soma and dendrites, and right through the axon, as depicted in Figure 2E. In the apical dendrite, the AP lasts considerably longer than in the soma and axon. More specifically, two AP spikes occur in the soma and axon during one spike in the apical dendrite for this somatic stimulation strength. During AP firing there are clear negative EP spikes near the soma and AIS, and the largest spikes are generated near the AIS, consistent with experimental observations in Bakkum et al. (2019). Figure 2D shows spatial solution snapshots of the membrane potential of the neuron and the surrounding EP near the soma and AIS during AP firing.

In our simulations, we consider three different types of stimulation.

I: Constant somatic stimulation

The first considered type of stimulation is a constant membrane current density, I_s (in μA/cm²), applied to the somatic part of the membrane. This membrane current density is included in the EMI model as described in (Equation 1). A positive I_s reduces the neuron's excitability, whereas a negative I_s increases the excitability of the neuron (see Figures 1C, 2C).

II: Constant extracellular stimulation

A second type of considered stimulation is an extracellular point source, representing extracellular stimulation through an electrode. This point source is represented by a sphere cut out of the mesh of the extracellular domain. On the surface of this sphere, a Neumann boundary condition of the form

is applied, where I_e (in μA) is the stimulation strength, A_e (in cm²) is the surface area of the sphere, and n_e is the outward pointing normal vector of the extracellular domain (i.e., pointing toward the center of the stimulation sphere). For a positive I_e, a negative EP is emanating from the stimulation sphere. The stimulation sphere is located in the vicinity of the AIS and the distance will be reported for each simulation. The EP surrounding the Purkinje neuron at rest as a result of a constant extracellular stimulation is illustrated for a few choices of I_e in Figure 3A. In addition, the EP at an AIS membrane point as a function of the electrode-to-AIS distance is provided in Figure 3B in a similar manner as in Figure 1B in Lee et al. (2024) and Figure 1C in Anastassiou et al. (2011).

III: Sinusoidal extracellular stimulation

The third considered stimulation type is a sinusoidal extracellular stimulation changing strength and sign with time during the simulation. This stimulation is implemented in the same manner as the constant extracellular stimulation, except that the boundary condition on the stimulation sphere is given by

where I_e (in μA) is the stimulation strength, A_e (in cm²) is the surface area of the stimulation sphere, f (in Hz) is the stimulation frequency, and t (in seconds) is time.

We solve the EMI system of equations numerically using the spatial and temporal operator splitting technique introduced in Jæger et al. (2021b) with a global time step of Δt = 0.01 ms. We apply the MFEM finite element library (Anderson et al., 2021; MFEM, 2026) with first-order elements and a mesh generated using Gmsh (Geuzaine and Remacle, 2009). For the non-linear membrane model dynamics, we use the first order Rush Larsen method (Rush and Larsen, 1978) with code generated using the Gotran software (Hake et al., 2020) and a time step of Δt = 0.001 ms for the Purkinje neuron and Δt = 0.01 ms for the pyramidal neuron. The numerical code is available on Zenodo (https://doi.org/10.5281/zenodo.18183929; Jæger, 2026).

As a first investigation into the effects of extracellular stimulation on neurons, we consider a non-firing cerebellar Purkinje neuron and a non-firing neocortical layer 5 pyramidal neuron located near a sinusoidal extracellular stimulation source with relatively weak stimulation amplitudes. The Purkinje neuron is made quiescent by a somatic stimulation current of I_s = 1 μA/cm², while no somatic stimulation is applied for the pyramidal neuron. The extracellular stimulation source is located 20μm from the AIS of the neurons. Figure 4 shows the membrane potential and the EP in an AIS membrane point of the two neurons. We observe that the extracellular stimulation yields an EP of up to about ± 0.5 mV and ± 1 mV for the stimulation amplitudes of I_e = 50 nA and I_e = 100 nA, respectively, at the AIS membrane. This results in a similarly sized alteration in the AIS membrane potential of the two cells. These subthreshold membrane potential fluctuations are present for different applied stimulation frequencies (1 Hz, 5 Hz, and 30 Hz), consistent with previous experimental findings (Anastassiou et al., 2011; Lee et al., 2024).

Next, we investigate how extracellular stimulation modulates the dynamics of firing neurons.

The next stage in our investigations is the study of ephaptic effects between two neurons. As a first example of such effects, we consider ephaptic synchronization of two firing cerebellar Purkinje neurons following up our preliminary analysis of this topic in Jæger and Tveito (2025b).

We now consider the case of ephaptic interactions between two neocortical layer 5 pyramidal neurons. In this case, we examine how the AP firing of a cell receiving substantial somatic stimulation modulates the AP firing of a nearby cell receiving weak somatic stimulation. Figure 10 shows such ephaptic interactions between two spiking pyramidal neurons. Like for the Purkinje neurons, we let the axon of the pyramidal neurons bend toward each other such that the minimum distance between the AIS of the two cells is 1 μm (see Figure 10A). The two cells are not connected by gap junctions or synapses. In Figure 10A, we show solution snapshots for a few time points during AP firing of the left cell (Cell 1). First, a negative EP is generated close to the AIS, and later toward the soma. Then, after a few milliseconds, a positive EP is present close to the soma and the AIS.

In Figure 10B, we illustrate the membrane potential at the AIS of Cell 1 for a number of different values of somatic stimulation (Is1) causing AP firing of different frequencies. As we vary the somatic stimulation of Cell 1, the somatic stimulation of Cell 2 is kept constant at Is2=-27μA/cm². The resulting Cell 2 AIS membrane potential is plotted in Figure 10C. In that simulation Is1=0. However, we want to investigate whether another choice of somatic Cell 1 stimulation, Is1, might alter the spike timing of Cell 2 even though the stimulation of that cell is kept the same. To that end, we define the spike distance for Cell 2 as illustrated in Figure 10C.

In Figure 10D, we report the spike distance of Cell 2 for different Cell 1 stimulation levels. We observe that as the stimulation strength of Cell 1 is increased and Cell 1 fires APs more frequently, the distance between spikes for Cell 2 increases by several milliseconds. We consider the default value of σ_e in addition to a reduced value of σ_e/5. The increase in Cell 2 spike distance as a result of rapid Cell 1 firing is more prominent for the reduced value of σ_e. Since the stimulation of Cell 2 is kept constant and the two cells are not connected by gap junctions or synapses, these spike time alterations must be caused by ephaptic interactions between the neurons.

As a final set of examples of ephaptic interactions between neighboring neurons, we consider the case of an AP being triggered in a neuron due to the EP generated from a nearby firing neuron. We consider the default value of σ_e = 3 mS/cm in addition to three cases of reduced conductivities. The setup used in the simulations is illustrated in Figure 11B.

For the pyramidal neurons, we let Cell 1 fire APs as a result of a somatic stimulation current of Is1=-50μA/cm². Cell 2 is given a weak somatic stimulation current of Is2=-25μA/cm², which is not strong enough to trigger AP firing in itself. In Figure 11A, we plot the membrane potential and EP in an AIS membrane point of each cell. We observe that for the default value of σ_e, APs are fired in Cell 1, but not in Cell 2. The EP generated during Cell 1 firing reaches a value of about −0.5 mV at the AIS of Cell 1 and about −0.2 mV at the AIS of Cell 2. The corresponding perturbation of the membrane potential, v = u_i−u_e, of Cell 2 is not large enough to cause AP firing. As the value of σ_e is reduced, the amplitude of the EP is increased, and we see clear perturbations of the membrane potential of Cell 2. Yet, even for σ_e reduced by a factor of 20, this perturbation is not large enough to trigger an AP. However, when σ_e is reduced by a factor of 50, the membrane potential perturbation of Cell 2 caused by the EP generated during Cell 1 firing is large enough to trigger AP firing at Cell 2 as well.

For the Purkinje neurons, Cell 1 does not receive any somatic stimulation current and fires spontaneous APs, while Cell 2 is given a positive somatic current of Is2=0.7μA/cm², which makes the cell quiescent on its own. However, if the extracellular conductivity is very low (default σ_e divided by a factor of 20 or 50), the negative EP generated from Cell 1 firing is enough to cause firing of Cell 2 as well (see Figure 11A).

A first step in examining how neurons can influence each other ephaptically is to understand how a single neuron responds to an external current source. This is also of independent interest, as it forms the basis for electrical stimulation techniques used in a range of experimental and clinical settings. The impact of an extracellular current source on a neuron has been investigated experimentally, see, e.g., Anastassiou et al. (2011) and Lee et al. (2024), and the spatial decay of the potential is known analytically from the solution of the electrostatic equation Holt and Koch (1999). If r denotes the distance in the extracellular space to a point current source, the resulting potential decays as 1/r, a pattern illustrated in Anastassiou et al. (2011) and Lee et al. (2024). The lower panel of Figure 3 follows this decay profile.

In Figure 4, a weak extracellular current source modulates the subthreshold membrane potential of a nearby neuron, and the resulting oscillation follows the sinusoidal frequency of the stimulus current. The amplitude of the membrane response increases in proportion to the strength of the applied current source. This type of subthreshold frequency-following, with an approximately linear dependence on stimulus amplitude, is consistent with experimental and modeling studies of neurons exposed to weak oscillatory electric fields (Anastassiou et al., 2011; Lee et al., 2024; Deans et al., 2007).

In Figure 5 we show how the spike frequency is affected by the strength of an external stimulation source and the distance to this source. For the Purkinje neuron, we first consider the case without somatic stimulation, i.e., I_s ≡ 0, see (Equation 1). In this setting, the number of spikes during a 500 ms simulation increases as the extracellular current becomes stronger and as the distance to the source decreases. However, for sufficiently strong stimulation at sufficiently close distances, the AP fails to repolarize properly, and the number of spikes is reduced to a single spike.

Next, we add a positive somatic membrane current to the Purkinje neuron, thereby reducing its excitability. This results in fewer spikes, while the dependence on the strength of the extracellular current source and the distance to the source remains similar to the case without added membrane current. Doubling the additional membrane current continues this trend and yields substantially fewer spikes.

The neocortical layer 5 pyramidal neuron is quiescent in the absence of membrane stimulation (i.e., I_s ≡ 0). In the lower part of Figure 5, we note that for this cell type the extracellular stimulus must be both close and strong in order to initiate spikes. When a sufficiently strong negative I_s is applied to the membrane, spikes are generated, and still more appear when the extracellular stimulus is close and strong.

In Figures 6, 7 we show how a sinusoidal extracellular current source affects the membrane potential and EP of the two neuronal models. For the Purkinje neuron, a weak low-frequency stimulus has limited impact on the membrane potential, but as the stimulation strength increases the timing of the AP firing is influenced in a phase-dependent manner by the sinusoidal drive. When the phase of the stimulation current produces a sufficiently positive perturbation in the EP, the membrane potential becomes silent and no spikes are generated, whereas the spiking frequency is increased in the opposite phase. These effects are also clear for a 5 Hz stimulation frequency. When the frequency is increased to 30 Hz or 100 Hz, low stimulation amplitudes again have limited impact on the spikes, whereas strong stimulation suppresses spiking during phases of strongly positive EP and activate spiking during phases of negative EP.

In Figure 7, the neocortical layer 5 pyramidal neuron is subjected to a weak somatic membrane stimulation in addition to the sinusoidal extracellular stimulation. We observe that the sinusoidal stimulation induces a combination of occasional spikes and subthreshold modulation rather than regular firing at the drive frequency. Nevertheless, spikes seem to preferentially be generated when the EP is negative corresponding to a negative phase of the stimulation current.

These findings are consistent with experimental studies showing that weak oscillatory electric fields modulate spike timing and entrainment in a cell-type-dependent manner (Lee et al., 2024) and that transcranial alternating current stimulation (tACS) can phase-lock single-neuron firing in vivo (Krause et al., 2019). They are also consistent with modeling work based on the Pinsky–Rinzel neuron model (Pinsky and Rinzel, 1994), which demonstrates changes in firing sensitivity and activity patterns under AC-induced electric fields (Yuan et al., 2025).

In both cardiac electrophysiology and neuroscience, the possibility that one cell can excite its neighbor without direct coupling via gap junctions or chemical synapses has been extensively discussed, see, e.g., Han et al. (2018); Blot and Barbour (2014); Anastassiou and Koch (2015); Bokil et al. (2001) for neural tissue and Veeraraghavan et al. (2014); Gourdie (2019); Jæger et al. (2023) for cardiac tissue. For both organs, this question goes to the core of how tissue organizes excitation waves and is therefore of considerable interest. Here, we assess ephaptic coupling for the two neurons under consideration.

In Figures 8, 9 we follow up the analysis presented in Jæger and Tveito (2025b), where we showed that two adjacent cerebellar Purkinje neurons synchronize using the modeling framework applied here. In Figure 8 we illustrate how the APs of two neighboring Purkinje neurons synchronize and how this process depends on the distance between the cells. When the cells are close (1 μm measured at the AIS), the synchronization is rapid and the two APs become completely synchronized. However, when the distance increases to 5 μm, the APs synchronize only partially and reach a stable state in which the spike firing of the two neurons is separated by about 1.2 ms. When the distance is increased to 10 μm, a similar behavior is observed, but the converged distance is about 1.8 ms. The emergence of near-synchronous firing between neighboring Purkinje cells in our simulations is consistent with in vivo recordings showing sub-millisecond, tightly synchronized activity of nearby Purkinje cells and population-wide spike alignment (Han et al., 2018; Sedaghat-Nejad et al., 2022), and with modeling studies demonstrating that decreasing inter-neuronal spacing or increasing packing density strengthens ephaptic interactions and promotes spike-time synchrony (Stacey et al., 2015; Goldwyn and Rinzel, 2016).

In Figure 9 we show how the speed of synchronization depends on the value of the extracellular conductivity, σ_e. Panel C demonstrates that synchronization is strongly affected by σ_e in the sense that reducing σ_e markedly increases the speed of synchronization. This is consistent with the observation in Tveito et al. (2017b) that the ephaptic current scales as O(1/σ_e). A similar sensitivity of ephaptic effects to the properties of the extracellular space has been reported in computational studies: in networks of ephaptically coupled neurons, stronger ephaptic potentials and reduced inter-neuronal spacing promote spike synchronization and microcluster formation (Stacey et al., 2015), in axon bundles, higher extracellular resistivity (equivalently, lower conductivity) together with increased fibre density enhance ephaptic coupling and facilitate near-synchronous spike volleys (Schmidt and R. Knösche, 2022), while earlier EMI simulations show that the ephaptic current scales approximately as O(1/σ_e) and thus becomes negligible for very large extracellular conductivities (Tveito et al., 2017b).

In Figure 10 we investigate ephaptic interactions between two neocortical layer 5 pyramidal neurons separated by 1 μm of extracellular space at the AIS. Here, Cell 2 is given a constant weak somatic stimulation current, while the somatic current of Cell 1 is stronger and varied between experiments. Panel A illustrates how an AP firing in Cell 1 perturbs the EP in the narrow cleft between the cells, potentially affecting the dynamics of Cell 2. Indeed, in Panel D, this effect is quantified by reporting the distance between two spikes of Cell 2 as a function of the somatic stimulation of Cell 1. The spike distance is observed to change by several milliseconds depending on the stimulation and resulting firing frequency of Cell 1. Furthermore, the effect increases when σ_e is reduced.

As mentioned above, a long-standing question in both cardiac electrophysiology and neuroscience is whether a cell can directly excite its neighbor through the extracellular space alone. We have seen above that one cell can influence another, but can this interaction trigger a full AP? A classical reference is Weingart and Maurer (1988), who examined impulse transfer between isolated myocyte pairs. When two single myocytes were gently pushed together in a side-to-side configuration, they observed no electrotonic interaction and no impulse transfer, and concluded that ephaptic transmission was unlikely. This experiment was revisited in Waaben et al. (2024), but with the cells aligned end-to-end rather than side-by-side. In this geometry, the high-density sodium channel regions at the intercalated disc of the myocytes are brought into proximity, and AP transmission was observed even when gap junctional coupling was absent. Thus, whether signal transfer occurs depends strongly on the cellular configuration.

In Figure 11 we explore this question for our two neuron types placed one μm apart. For the EP generated by one cell to be strong enough to elicit a full AP in its neighbor, the extracellular conductivity σ_e must be sufficiently small. This is consistent with the ephaptic current scaling as O(1/σ_e). Whether such small σ_e values occur locally in neural tissue remains uncertain, although geometric confinement and narrow extracellular clefts could, in principle, create domains with reduced effective conductivity. However, in Figure 11, we find that for the EP generated by one cell to directly trigger an AP in a neighboring neuron separated by one μm, the extracellular conductivity, σ_e, must be reduced by roughly a factor of twenty to fifty relative to the standard value of 3 mS/cm. We repeated the simulations with the distance between the neurons reduced by one half, but this did not change the conclusion regarding the requirements for direct excitation. This places the direct triggering regime outside the range of physiological conductivities. In healthy neural tissue, bulk extracellular conductivity is typically in the range of 3–6 mS/cm (Logothetis et al., 2007).

A notable exception to these bulk values is the Pinceau structure surrounding the Purkinje cell AIS, which forms a locally high-resistivity compartment and supports strong ephaptic inhibition (Blot and Barbour, 2014). However, this structure is anatomically and functionally specialized, and its ephaptic effect is inhibitory rather than excitatory. Thus, while Figure 11 shows that direct ephaptic triggering is theoretically possible in the extreme limit of conductivity, it is unlikely to operate as a mechanism for excitatory propagation in typical cortical or cerebellar tissue. Instead, ephaptic effects are more plausibly involved in modulating spike timing and supporting synchronization, as observed at physiological conductivities in Figures 9, 10.

Several numerical implementations of the EMI system (1) have been proposed. In Tveito et al. (2017b), a straightforward finite difference scheme was used. The complete domain was represented by a finite difference grid in which intra- and extracellular nodes overlap at the membrane, allowing all unknowns to be fully coupled. A finite element approach was presented in Tveito et al. (2017a) and re-used in Buccino et al. (2019). Furthermore, a boundary element approach was used in de Souza et al. (2024), and a cut finite element method was used in Berre et al. (2024).

The implementation used in the present project is described in detail in Jæger et al. (2021b). The underlying idea is to use operator splitting to formulate the problem in terms of standard components in scientific computing: Poisson/Laplace problems in the intracellular and extracellular domains, coupled to ordinary differential equations at the membrane nodes. This decomposition enables straightforward use of well-tested solvers from standard numerical libraries. In this work we use the MFEM finite element library, but any standard finite element library will work well, since the spatial subproblems are standard Poisson/Laplace equations; the key point is to reduce the coupled EMI system to these standard components, as described in Jæger et al. (2021b). The applied numerical code is available on Zenodo (Jæger, 2026).

For both cardiomyocytes, see, e.g., Weingart and Maurer (1988); Waaben et al. (2024); Veeraraghavan et al. (2014); Gourdie (2019), and neurons, see, e.g., Anastassiou et al. (2011); Anastassiou and Koch (2015); Han et al. (2018), there has been a lengthy debate on whether ephaptic coupling is of physiological significance. One reason for this debate is that ephaptic effects are expected to be strongest in strictly confined extracellular spaces, where direct measurements of voltages are technically challenging. Computational models have therefore been commonly applied (Sperelakis and Mann, 1977; Schmidt and R. Knösche, 2022; Jæger et al., 2023), but these studies have often lacked accompanying experimental validation.

For the present study, partial validation of the overall EMI framework could be obtained by comparing stimulation–response curves to experiment. For example, the spike counts presented in Figure 5 for both the Purkinje neuron and the neocortical layer 5 pyramidal neuron are strongly modulated by the distance to the external current source and by the stimulation strength. Such dependencies are directly measurable in vitro (or ex vivo) using controlled extracellular stimulation, and could be compared to model predictions in terms of spike count, spike timing shifts, and threshold currents as functions of stimulation amplitude and distance.

Although the present computational scheme based on the EMI model allows for a more detailed representation of neurons than commonly used cable-based models, there are several important limitations. First, the 3D neuron geometries were generated in Gmsh as idealized morphologies and were not constrained by anatomical reconstructions. Second, the membrane models contain a long list of parameters that are guided by experimental measurements, but these parameters are associated with substantial uncertainty and non-uniqueness, in the sense that quite different parameter sets may reproduce similar activity patterns (Marder and Taylor, 2011; Prinz et al., 2004). Third, we have assumed that ionic concentrations are constant in space and time within both the intracellular and extracellular domains. These concentrations could have been allowed to vary by, e.g., applying the KNP-EMI model, which combines the EMI framework with an electroneutral electrodiffusion description and allows ion concentrations and reversal potentials to evolve dynamically (Ellingsrud et al., 2021, 2020). Moreover, at nanometer scales near membranes and in narrow extracellular clefts, the homogenized description used in KNP-EMI breaks down, and electrodiffusive effects, Debye layers and steep local concentration gradients become important. Such effects can be captured only by nanoscale Poisson-Nernst-Planck models, see, e.g., Pods (2017); Pods et al. (2013); Jæger and Tveito (2025a), and are not included in the EMI model formulation applied here.