A Novel Neural Computational Model of Generalized Periodic Discharges in Acute Hepatic Encephalopathy

A1). Overview of the Liley model

The Liley model comprises two neural populations: an excitatory population (represented by E) and an inhibitory population (represented by I), which are functionally distinct but synaptically coupled. The model includes two inter-connections (weighted by N_ie and N_ei) and two intra-connections (weighted by N_ee and N_ii). External input from thalamus to E and I are presented by p_ee and p_ei. Fig. 1 illustrates the topological structure.

The dynamics of each population are described by two state variables, the soma membrane potential and the post-synaptic potential. The soma membrane potential V is modeled according to conductance-based rules [24], specifically,

where

Here, V^r is the resting membrane potential and V^eq represents the reversal potential. On the other hand, post-synaptic potential I is modeled according to the convolution-based rules, specifically,

where

Note that m(x, t) and h(t) denote the firing rate and impulse response of a population where t represents time and x corresponds to the spatial position of the population in the brain. In (4), S(∙) is a sigmoid function, Q_max is the maximum population firing rate, σ is the steepness of the sigmoidal transformation, and V₀ is firing threshold. In (5), Γ represents the peak amplitude of the post-synaptic potential and γ characterizes the exponential decay time scale of the post-synaptic potential. For convenience in calculating the convolution (⊗ in (3), further taking the Laplace transformation on it, then Eq. (3) can be re-expressed by a non-linear partial differential equation as

Details of the transformations of electrical activities in the Liley model are shown in Fig.2. Blue lines represent intra-connected loops and the green dotted lines represent inter-connected loops. In each loop, the post-synaptic potential I can be obtained from the soma membrane potential V by Eq. (3), which is completed by two computational blocks (see sub-figure in Fig. 2): a potential-to-rate block (represented by “S”) and a rate-to-potential block (represented by “h”). Specifically, the soma membrane potentials V is transformed into firing rate m by Eq. (4) in the first block, and the firing rate m is converted into the post-synaptic potential I by Eq. (5) in the second block. After that, the post-synaptic potential I is further transferred into the soma membrane potential V according to Eqs. (1)-(2).

A2). Neuropathological alterations of AHE

In acute liver failure, there is a systemic build-up of toxic substances from the gut, including ammonia, that are normally metabolized and excreted by the liver [3]. These toxins travel via the bloodstream to the brain, ultimately leading to development of AHE. Among numerous cellular and neurophysiological mechanisms underlying AHE, we we focus here on three for which there is clear evidence for a role in development of AHE. These are (1) impaired synaptic transmission secondary to reduced cellular metabolism; (2) increased neuronal excitability; and (3) enhanced postsynaptic inhibition.

Brain concentrations of ATP are significantly reduced in experimental animal models of AHE [25][26]. This is thought to result from hyperammonia-induced mitochondrial dysfunction, systemic ischemia and hypoxia and increased ATP consumption [25][27]. The latter process may be secondary to hyperammonia-induced glutamate N-methyl-D-aspartate (NMDA) receptor hyperactivation and concomitant compensatory excessive activation of the sodium/potassium ATPase pump to maintain cellular sodium levels within the normal range [25][27]. Given that synaptic neurotransmission is a highly metabolically expensive process, accounting for 30% of the brain’s adenosine triphosphate (ATP) usage [28], it is likely that ATP depletion observed in AHE is associated with impaired synaptic signaling. Supporting this claim, studies have shown that the reduction in brain ATP levels coincides with the onset of AHE-associated behavioral impairment and EEG abnormalities [25][26].

Increased neuronal excitability is another key mechanism underlying AHE. Increased neuronal excitability in AHE appears to arise from two processes: hyperactivation of glutamate NMDA receptors and anoxic long-term potentiation. Ammonia levels are significantly increased in AHE patients [29], and elevated ammonia levels in acute liver failure (ALF) animal models lead to overactivation of NMDA glutamate receptors ([30]). Enhanced NMDA-mediated neuronal excitability is also supported by studies which demonstrate decreased expression of astrocytic glutamate transporters (EAAT-2) [31] and astrocytic glycine transporters (GLYT-1) [14], corresponding to increased levels of synaptic glutamate and glycine, agonists and co-agonists of NMDA receptors respectively, in ALF rat models and AHE patients. The resultant overactivation of glutamate NDMA receptors is followed by excessive formation of nitric oxide and cyclic guanosine monophosphate (cGMP), and ultimately increased neuronal excitability [30]. In addition, anoxic long-term potentiation, the phenomenon of increased excitatory postsynaptic potential (EPSP) amplitudes following exposure to ischemia or anoxia, may also contribute to increased neuronal excitability in AHE. Indeed, AHE is often accompanied by a variety of other medical processes which may lead to decreased brain oxygenation, such as gastrointestinal bleeding, sepsis, or the effects of cytokines and other toxins released from necrotic liver tissue [32][33].

Finally, a growing body of evidence suggests that enhanced gamma-aminobutyric acid (GABA)- and glycine-mediated postsynaptic inhibition may also be central to the pathophysiology of AHE. Levels of pregnenolone-derived neurosteroids (e.g. allopregnanolone), potent selective positive allosteric modulators of the GABA_A receptor, are significantly increased in experimental ALF models and in brains of hepatic coma patients [34]. Additionally, a study in a rabbit model of ALF has shown greatly increased levels of GABA-like activity in peripheral blood plasma just before the onset of AHE [35]. As with ammonia, this finding may be partially due to impaired hepatic extraction of gut derived GABA from portal venous blood [35]. The same study also showed that the blood-brain barrier becomes abnormally permeable to an isomer of GABA before the onset of AHE, and that hepatic coma is associated with an increase in the brain density of GABA receptors [35], thereby increasing the sensitivity of the brain to GABA-ergic neural inhibition. Another study showed that treatment with flumazenil, a GABA-benzodiazepine receptor antagonist, appears to result in clinical and electroencephalographic improvements in AHE patients [36]. Further, increased levels of synaptic glycine, an inhibitory neurotransmitter observed in AHE, may lead to over-activation of glycine receptors, which together with the hyperactivation of GABA receptors result in post-synaptic inhibition through chloride-channel opening.

A3). Neural computational model of AHE

To begin building our neural computational model of AHE (AHE-CM), we mathematically augment the Liley model to capture the three pathophysiologic mechanisms described above for AHE. We consider only temporal dynamics; we do not model spatio-temporal dynamics in this work.

First, the postsynaptic potential amplitude Γ in Eq. (6) is taken as a variable satisfying the slow dynamics of itself coupled to the presynaptic firing rate S(∙) as in the work of Bojak et.al. [37], which can be expressed by

Here, $p_k^{dep}$ is the depletion constant reflecting the impaired synaptic transmission, ${\Gamma }_k^{rest}$ is the resting value of the maximum amplitude of the postsynaptic potential (PSP), and ${\tau }_k^{rec}$ is the recovery time for activity dependent synaptic depression, subscript k ϵ {e, i} indicate excitatory and inhibitory populations, respectively. Note that this mechanism applies to both excitatory and inhibitory populations.

Next, to model the effect of increased neuronal excitability, we modify the Liley model to amplify the resting value of maximum amplitude Γ^rest of the excitatory postsynaptic potential (EPSP) in Eq. (7). This is accomplished by setting,

where ${\Gamma }_e^{equ}$ is the equilibrium voltage of the EPSP and F_am is the amplification factor. This mechanism applies only to the excitatory population.

We last model enhanced postsynaptic inhibition, caused by the increased GABAergic neural transmission. To do this, we prolong the duration of the inhibitory postsynaptic potential (IPSP) by altering the decay rate constant γ_i in Eq. (6) according to [38], so that the inhibitory post-synaptic potential is given by

where

Here, ${\gamma }_i^0$ is the baseline synaptic time constant, ϵ_i is the control parameter relating to the decay time of the IPSP. This mechanism applies only to the inhibitory population.

On the basis of the dynamics of each population and detailed transformation of electrical activities in Liley model, the proposed neural computational model of AHE (AHE-CM) can be formulated by combining the mathematical expressions (7)-(9) of three mechanisms. Eq. (11) presents the mathematical model of AHE-CM. The Euler-Maruyama method with a time step of 0.1ms is applied to solve the Eq. (11) and the soma membrane potential of the excitatory population (i.e. V_e) is taken to represent the EEG signal. We denote by Θ the set of parameters in AHE-CM. A detailed explanation of the biological interpretation of these parameters and their default values, which we adopt, can be found in [38].

Modeling impaired synaptic transmission

To characterize impaired synaptic transmission in AHE, we follow the idea in the bursting Liley model of including synaptic depression [37], which means that the postsynaptic peak amplitudes Γ_e and Γ_i decrease as a function of presynaptic firing rates and recover with time constants ${\tau }_e^{rec}$ and ${\tau }_i^{rec}$ (i.e. Eq. (7)). We set baseline values of ${\tau }_e^{rec}$ to 500ms and ${\tau }_i^{rec}$ to 250ms, in agreement with the physiological range of 250ms to 1000ms [39]. According to the work of Ruijter et al.[23], the impaired synaptic transmission will induce ${\tau }_e^{rec}$ and ${\tau }_i^{rec}$ that are dramatically increased. Therefore, ${\tau }_e^{rec}$ and ${\tau }_i^{rec}$ are varied from their baseline to their maximum value. When varying these parameters in our experiments, because excitatory synapses are believed to recover more slowly than inhibitory synapses [40], we allow a larger range for ${\tau }_e^{rec}$ [500 20000]ms than for ${\tau }_i^{rec}$ [250 2500]ms, as in [23].

Modeling elevated neuronal excitability and enhanced postsynaptic inhibition

To model the effect of elevated neuronal excitability, we increased the resting value of maximum EPSP (${\Gamma }_e^{rest}$) through the parameter “F_am”. Note that, as shown in Eq. (7), Γ_e is reduced because of the impaired synaptic transmission and may recover with time constant ${\tau }_e^{rec}$ to its baseline value. While through the parameter F_am, it may increase above its baseline value, reflecting increased neuronal excitability. We set the magnitude order of F_am to the same values used in the work of Rujiter et.al [23] in modeling HIE.

To model enhanced gamma-aminobutyric acid and glycine-mediated postsynaptic inhibition, we followed the approach of [38], which modeled effects of anesthesia with GABA_A agonists (e.g. propofol). The baseline value of the synaptic rate constant (${\gamma }_i^0$) is 0.065ms⁻¹, which agrees with the physiological range of 0.01ms⁻¹ to 0.5ms⁻¹ [23]. The magnitude order of ϵ_i is chosen according to [41]. The postsynaptic potential decay rate γ_i can be calculated by Equation (10).

B1). Performance verification of three added mechanisms

Here we study model behaviors in response to changes of key parameters, and illustrate the effects of three extra mechanisms added to the Liley model. The four parameters ${\tau }_e^{rec}$, ${\tau }_i^{rec}$, F_am, γ_i, which are employed in the three mechanisms are considered the key parameters in the proposed AHE-CM. All other parameters in the model are set to their nominal values in Liley et al. [46], summarized in Table 2. The external input p_ee is generated randomly based following prior literature, by sampling from a Gaussian distribution with mean 3460s⁻¹ and standard deviation 1000s⁻¹; p_ei is constant (5070s⁻¹) [23]. We first focus on studying qualitatively what kind of model outputs are produced by different combination of key parameters. In all simulations, EEG signals of 50s duration were generated; the first 10s were ignored to exclude transient effects. The Matlab code used for the model simulations is provided in https://github.com/jill-Song/AHE-code.

The 7 types of simulated signals of AHE-CM are categorized according to the mean value of the inter-TPWs intervals (M-ITIs) in the simulated EEG segment (shown in Table 3). In Fig. 5 (a), we define a contrast between the 7 types of model outputs and their labels. The model output with label 1 represents a “low voltage” pattern. Labels 2 to 7 exhibit periodic discharges of increasing frequencies. Fig. 5 (b) illustrates the results under 9 different parameter settings, where each pair (F_am, γ_i) is taken from (3, 1.5, 0.8) × (0.065, 0.032, 0.014). In each case, ${\tau }_e^{rec}$ varies from 500 to 20000 with step-size 100, and ${\tau }_i^{rec}$ varies from 250 to 2500 with step-size 50. That is, each graph depicts 195 × 45 labels of model outputs. From Fig.5 (b), it can be seen that higher frequency periodic discharges occur for larger values of (F_am, γ_i), and smaller frequency discharges occur in the plots in the right panels with lower values of these parameters. On the other hand, in a graph with fixed F_am and γ_i, the discharges increase in frequency (label numbers increase) as ${\tau }_e^{rec}$ decreases and as ${\tau }_i^{rec}$ increases. These trends emerges very clearly in “graph (3)” where almost all patterns/labels occur.

We next examine one special case (as shown in Fig. 6 (d)) in detail to illustrate the effects of adding three extra mechanisms added to the baseline Liley model. Fig. 6 illustrates the outputs of four different models, each with one of the additional mechanism added in. In Fig. 6 (a), the output of the baseline Liley model is shown, where all parameters are set to their nominal values (see Table 2). Fig. 6 (b) illustrates the output of Liley model with the first mechanism added, where the postsynaptic potential amplitude Γ is determined by Eq. (7) with ${\Gamma }_k^{rest}=0.71(k\in \{e,i\})$, ${\tau }_e^{rec}=6000$, ${\tau }_i^{rec}=1000$. We see that periodic waves begin to appear. The periodicity becomes more evident in Fig. 6(c), as the second mechanism is added, that is, when Γ is determined by both (7) and (8); in this illustration F_am = 0.8. Finally, the output of AHE-CM is obtained by adding all of three mechanisms, shown in Fig. 6 (d). Here, except for the estimated value of Γ as mentioned above, the decay rate constant γ_i is further calculated by (10). Interestingly, although reproducing the morphology of TPWs was not a goal of our modeling work, the close up in Fig. 6(d) shows that the output of the AHE-CM does provide an approximation of the morphology of TPWs. Each of the three red dots in close up of the waveform corresponds to one of the three phases of a TPW.

B2). Performance verification of AHE-CM combining with PF-POM

Next, we verify the ability of the proposed PF-POM method to infer the parameters of AHE-CM. We do this in simulated EEGs (generated by the AHE-CM model), because for these we know the true values for the underlying parameters. These true values can be served as a ground truth.

For simplicity, we consider only two key parameters ${\tau }_e^{rec}$, ${\tau }_i^{rec}$. We set F_am = 0.8, γ_i = 0.032, and all other parameters are set to their be nominal values in Table 2. According to the mentioned previously, the ranges of ${\tau }_e^{rec}$ and ${\tau }_i^{rec}$ are suggested to be [500,20000]ms and [250, 2500]ms. We first divide the parameter space into 5 sub-spaces by segmenting the interval of ${\tau }_e^{rec}$ and keeping the interval of ${\tau }_i^{rec}$ whole (see Table 4). Given an AHE-EEG signal S, we calculate the sequence of inter-TPW intervals $\mathcal{I}(\mathcal{S})$ as well as its mean, denoted by mean($\mathcal{I}(\mathcal{S})$). The sub-space related to S is then selected according to the value of mean($\mathcal{I}(\mathcal{S})$) shown as in Table 4. With the selected sub-space, the total number of particles is set to be N = 35 in our simulations.

In our PF-POM, the feature F₁ is defined to be the frequency of TPWs in AHE-EEGs. To illustrate this feature, Fig. 7 shows the frequency histogram of $\mathcal{I}$ with respect to two signals (S_fast and S_slow) with different frequencies (fast and slow) of TPWs. In Fig. 7 (b), the frequency histogram of $\mathcal{I}(S_{fast})$ is left-skewed normal distribution toward lower inter-TPW intervals, reflecting the higher average frequency of simulated TPWs signal S_fast. By contrast, the histogram of $\mathcal{I}(S_{slow})$ are more uniformly distributed and are concentrated at higher inter-TPW intervals, reflecting the slower average rate of TPWs in S_slow.

We next verify the ability of PF-POM to estimate model parameters from simulated EEG data. Given parameters values $({\tau }_e^{rec},{\tau }_i^{rec})=(1100,600)$, the model output (“simulated EEG”) is computed according to Eqs. (11). Applying the PF-POM yields estimates $({\tau }_e^{rec}{)}^{\ast }$ and $({\tau }_i^{rec}{)}^{\ast }$. Fig. 8 illustrates the results for 50 trials, where the mean and standard variation of them are shown. To further verify the performance of PF-POM, we compare model output with the true underlying parameters to the model output using the estimated parameter values. In Fig. 9, the red curve shows the original simulated EEG signal (labeled as “ground truth”), and the green curve shows the model output with the optimal values determined by PF-POM (here, one of the results in 50 trials is randomly selected to create the green curve and labeled as “simulated”). The comparison illustrates that effectiveness of the proposed PF-POM.

B3). Relation between features of AHE-EEGs and the mechanistic model of AHE

In this section, we focus on the relation between the frequency of TPWs in real AHE-EEGs and the mechanisms underlying AHE in our computational model. We have seen in Fig. 5 that variation in TPWs’ rate depends mainly on ${\tau }_e^{rec}$ and ${\tau }_i^{rec}$. Especially, when F_am = 0.8, γ_i = 0.065, the much more clear variation in all TPW patterns can be observed. Therefore in what follows, we use PF-POM to optimize only ${\tau }_e^{rec}$ and ${\tau }_i^{rec}$ (F_am = 0.8, γ_i = 0.065), in the AHE-CM, and leave the values of the other variables at their nominal values.

In our analysis of AHE-EEG data from patients, each signal is divided into 20s EEG segments. Then given each AHE-EEG segment, we use PF-POM to estimate the values ${\tau }_e^{rec}$, ${\tau }_i^{rec}$ of an AHE-CM model that best match the data. We then use the estimated AHE-CM model generate a simulated EEG segment (i.e., output of AHE-CM). Table 5 shows the number of TPWs in AHE-EEGs (TA) and simulated EEGs (TS), the mean of ITIs in AHE-EEGs (MA) and simulated EEGs (MS), as well as the corresponding “errors” / differences. Comparing “No. of TPWs” in 30min-long AHE-EEGs and simulated EEGs, the maximum error is 0.0734 and the minimum error is 0.0322, obtained by $\frac{\mid TA-TS\mid }{TA}$. Moreover, the errors between “mean of ITIs” in AHE-EEGs and simulated EEGs, |MA – MS|, fall inside the range [0.01, 0.06]. We thus conclude that the simulated EEGs resulting from our proposed approach “AHE-CM + PF-POM” have similar properties to AHE-EEGs with respect to periodicity of TPWs. Fig. 10 illustrates a 20s AHE-EEG segment from “patient 1” and the simulated EEG segment from the corresponding AHE-CM model after parameter optimization. The pink columns correspond to the annotated TPWs. We observe that 30 TPWs are included in the AHE-EEG segment, whose mean of ITI is 0.6442 (Fig. 10(a)). The model with estimated parameters produces 29 TPWs with mean ITI 0.6679 (Fig. 10(b)). These results suggest that optimizing ${\tau }_e^{rec}$ and ${\tau }_i^{rec}$ indeed is able to qualitatively match the periodic discharge behavior of real AHE EEGs.

We next focus on characterizing the relation between the rate of TPWs in AHE-EEGs and the mechanisms underlying AHE. From the above experimental results, we select 9 out of 90 segments for each patient (shown in the Supplementary material), and the corresponding mean IEI values and optimal parameter settings (${\tau }_e^{rec}$, ${\tau }_i^{rec}$) estimated for each segment. In Fig. 11, 63 points (9 × 7) are marked whose position corresponds to the optimal parameters determined by PF-POM and whose value is the mean of its ITIs (M-ITIs). We observe that EEG segments having comparatively small values of M-ITIs (i.e. high frequency of TPWs) locate at the bottom right corner of the graph (corresponding to small ${\tau }_e^{rec}$ and large ${\tau }_i^{rec}$), and EEG segments are located at the upper-left of the graph having comparatively high values of M-ITIs (corresponding to small ${\tau }_e^{rec}$ and large ${\tau }_i^{rec}$).

These results are plausible from a physiological point of view. Here, parameters ${\tau }_e^{rec}$ and ${\tau }_i^{rec}$ represent the recovery times for excitatory and inhibitory activities respectively. Therefore, the less the value of ${\tau }_e^{rec}$ is, the more quickly recovery of excitatory activity will be, thus the faster the excitatory activity cam resume - leading to higher frequency discharges. Similarly, the excitatory activity can also recover more quickly if the value of ${\tau }_i^{rec}$ is larger. On the basis of the above analysis, we hypothesize that the discharge rate of TPWs in real AHE-EEGs is related to recovery times of excitatory and inhibitory neural populations.