Introduction

Normal aging often leads to cognitive impairment, as evidenced by reduced performance on learning and memory tasks in both humans (Albert, 1993; Salthouse et al., 2003; Fisk and Sharp, 2004; Rhodes, 2004; Sorel and Pennequin, 2008) and non-human primates (Moore et al., 2006; Shamy et al., 2011; Moore et al., 2017; Comrie et al., 2018; Chang et al., 2022; Moore et al., 2023). In the rhesus monkey, age-related working memory decline is accompanied by sublethal structural and functional changes in vascular elements, individual pyramidal neurons, glial cells, and white matter pathways (Hof and Morrison, 2004; reviews: Luebke et al., 2010 and Peters and Kemper, 2012; Morrison and Baxter, 2012; Luebke et al., 2015; Shobin et al., 2017; Dimovasili et al., 2022). It is well documented that neurons do not die during normal aging but rather undergo a number of morphological and physiological alterations, particularly in the monkey dorsolateral prefrontal cortex (dlPFC), the critical cortical circuit for working memory. For example, during normal aging Layer 3 pyramidal neurons in rhesus dlPFC exhibit a significant loss of dendritic spines and synapses (Chang et al., 2005; Peters et al., 2008), and electrophysiological changes observed both in vitro (Chang et al., 2005; Ibañez et al., 2020) and in vivo during a spatial working memory task (Wang et al., 2011). Perhaps most strikingly, extensive myelin dystrophy during normal aging has been observed in both gray and white matter (Peters et al., 2001; review: Peters, 2007; Bowley et al., 2010), including in monkey dlPFC (Peters and Sethares, 2002, 2003; review: Peters, 2009). Ultrastructural studies reveal that 3-6% of myelin sheaths in dlPFC exhibit age-related alterations including splitting of the major dense line of the myelin sheath, balloons, and redundant myelin (Figure 1; Peters and Sethares, 2002). Remyelination has also been observed across the adult lifespan in monkey dlPFC. Electron microscopy (EM) revealed a 90% increase in the number of paranodal profiles in aged versus young monkeys, indicating higher numbers of internodal myelin sheaths with aging. Aged subjects also had a significant proportion of abnormally short and thin myelin sheaths (Peters and Sethares, 2003). The hypothesized mechanism to explain these findings (review: Peters, 2009) is that myelin degradation begins as oligodendrocytes degenerate due to oxidative stress, and that axons accumulate dense inclusions in spaces between the lamellae of their associated myelin sheaths. As oligodendrocytes die, the associated sheaths detach from the axolemma, leaving bare axonal sections (complete demyelination). Subsequently, new mature oligodendrocytes remyelinate the bare segments, but with shorter and thinner sheaths. It is believed that the altered sheaths lead to a slowdown of signal propagation that contributes to cognitive impairment with aging. Several of the changes that pyramidal neurons undergo with aging correlate with observed cognitive impairment (review: Luebke et al., 2010; Peters and Kemper, 2012; Shobin et al., 2017; Moore et al., 2023), including myelin dystrophies and remyelination (Peters and Sethares, 2002, 2003; Dimovasili et al., 2022). However, which changes are the key determinants of age-related cognitive decline has not yet been firmly established (Konar et al., 2016; Motley et al., 2018; Cleeland et al., 2019). This is in part due to the difficulty of isolating individual neuronal features (e.g., firing rate, synapses, myelin) empirically, while controlling for concomitant changes in others. Thus, computational models become essential to predict how age-related changes in individual neurons affect cognitive impairment.

Figure 1.

Numerous modeling studies have examined connections between axon parameters and action potential (AP) conduction velocity (Rushton, 1951; Goldman and Albus, 1968; Brill et al., 1977; Moore et al., 1978; Waxman, 1980; Chomiak and Hu, 2009). Others have sought the most appropriate way to model the axon (e.g. Blight, 1985; Richardson et al., 2000; McIntyre et al., 2002; Gow and Devaux, 2008; Dekker et al., 2014). Demyelination has been modeled frequently in the context of disease (review: Coggan et al., 2015), showing that loss of myelin leads to slower AP conduction velocity (CV) and sometimes to AP failure. Remyelination of axons with shorter, thinner internodes decreases CV in models of both large- and small-diameter axons (Lasiene et al., 2008; Powers et al., 2012; Scurfield and Latimer, 2018). However, to our knowledge, the effects on AP propagation of widespread demyelination and remyelination in a broad population of axons has not been explored systematically. There have also been studies of how CV and myelin alterations impact network synchrony in phase oscillator models (Karimian et al., 2019; Pajevic et al., 2014; Noori et al., 2020; Pajevic et al., 2023), but not in spiking neural networks that can explain the mechanisms underlying working memory. Our recent network model found that age-related changes in firing rates and synapse numbers in individual neurons can lead to working memory impairment (Ibañez et al., 2020), but did not consider myelin dystrophy.

This study unites models at two different scales – multicompartment models of dlPFC pyramidal neurons and a spiking network model simulating spatial working memory – to investigate how age-related myelin degradation (represented by demyelination) and remyelination affect signal transmission and working memory precision. Higher degrees of demyelination led to slower propagation and eventual failure of APs along the axons, resulting in working memory impairment in the network model. Sufficient remyelination led to recovery of signal transmission and working memory performance. Overall, our study indicates that observed levels of myelin changes, if uncompensated by other factors, would lead to substantial working memory impairment with aging.

Results

Progressive demyelination causes CV slowing and AP failures in model neurons

To simulate myelin alterations in individual neurons, we adapted our multicompartment model tuned to data from rhesus monkey dlPFC (Rumbell et al., 2016) by attaching an axon model that captured nodes and detailed myelinated segments (Gow and Devaux, 2008; Scurfield and Latimer, 2018). Myelinated segments included an internode with adjacent juxtaparanodes and paranodes, and tight junctions between the innermost myelin lamella and axolemma (see Methods; Figure 2A-B). We applied demyelination and remyelination perturbations to a cohort of 50 young (control) neuron models, with axonal parameters varying within biologically plausible ranges. The CV in control models varied across the cohort and in response to myelin alterations (Figure 2C; Supplementary Figure 1). CV changes were more sensitive to variations along some dimensions of the parameter space than to others (e.g. myelinated segment length versus axon diameter), explored further below.

Figure 2.

AP propagation was increasingly impaired with increasing levels of demyelination across the control model cohort (Figure 3A): CV became slower, eventually leading to AP failure in several cases. Summarizing these results as a function of segments affected and lamellae removed made the overall trends clearer (Figure 3B). For example, irrespective of the extent of demyelination along an axon, when only 25% of lamellae were removed, CV change was negligible. However, when 100% of lamellae were removed, CV slowed drastically – by 38 ± 10 % even when just 25% of the segments were demyelinated. Under similar circumstances, when 75% of the segments were demyelinated, mean CV change was as much as -72 ± 8 %. Indeed, several axons under this condition showed AP failures, as evident from heat maps in Figure 3A. The percentage of high failure trials, in which at least 75% of control APs were lost, increased dramatically as demyelination increased (Figure 3C).

Figure 3.

Remyelination causes partial recovery from CV slowing and AP failures

We next examined how remyelination, occurring after demyelination, affected axonal AP propagation in the neuron models (Figure 4). We first assumed that affected segments had been previously completely demyelinated, i.e. losing all their lamellae (Figure 4A-C). Most remyelinated models showed CV recovery from 0-100%, except for a few cases in which CV increased relative to the unperturbed models (Supplementary Figure 2; see Discussion). Regardless of the extent of initial demyelination along an axon, CV recovery increased as both the remyelination and the lamellae restoration percentages increased. When all demyelinated segments were subsequently remyelinated – so that none of the perturbed segments were bare (unmyelinated) – with sufficient myelin lamellae, conduction recovery was similar to that in the unperturbed cohort (Figure 4B). A milder dependence of CV recovery on the initial fraction of demyelination also emerged. When all demyelinated segments were remyelinated, there was a positive relationship between the initial demyelination rate and the CV recovery: CV recovered more when 75% of the segments were demyelinated (Figure 4B, black lines). This finding is plausibly explained by observations of Scurfield and Latimer (2018), in which remyelinated axons with a larger number of transitions between long (unperturbed) and short (remyelinated) segments had slower CV. We saw the same relationship in our simulations (Supplementary Figure 3). When remyelination was less than 100%, there was a negative relationship between the initial amount of demyelination and CV recovery, explained by the number of bare segments that were not remyelinated. Results for the percentage of high failure trials (Figure 4C) were consistent with those for CV recovery. While the percentage of high failure trials decreased across the cohort as remyelination increased, certain models were more prone to AP failures than others, motivating the statistical analysis below. Consequently, the high failure trials percentage was more variable (wider error bars) than the CV recovery.

Figure 4.

We also simulated remyelination after partial demyelination, where affected segments initially lost only half their lamellae (Figure 4D). Overall trends were similar to those for the complete demyelination case (Figure 4E-F; Supplementary Figure 4). As expected, these simulations showed more CV recovery than when starting with complete demyelination, with less separation between the curves for different levels of remyelination. The variability in CV recovery across different remyelination conditions was similar to the variability across the 50-neuron cohort, especially when remyelinated segments had less than 50% of lamellae restored (i.e., less than the partially demyelinated segments). There were also fewer high failure trials under partial demyelination conditions, relative to the corresponding complete demyelination cases (Figure 4C vs. 4F).

Several parameters contribute to CV changes after myelin perturbation

Taken together, our multi-compartment neuron model simulations showed that both demyelination and remyelination can cause substantial CV delays and AP failures within individual axons. However, responses to the same simulated dystrophy sometimes varied widely across the cohort. To identify the key parameters that contribute to CV changes, we employed lasso regression separately on the demyelination and remyelination results. For CV changes in response to demyelination (Figure 5A), only 5 out of 12 parameters had non-zero coefficients in the lasso. Among these, the weight for myelinated segment length had the largest magnitude (and was negative). This is intuitive: models with longer myelinated segments show more CV slowing in response to a given demyelination perturbation. Scale factors for leak and sodium conductance, axoplasm resistance, and tight junction resistance also affected CV change substantially in the demyelination condition. For CV recovery after remyelination (Figure 5B), the lasso selected a model with 10 of the 12 parameters (only myelin length and axoplasm resistance were omitted). To test the effectiveness of lasso, we utilized these reduced models to predict how CV would change after demyelination and remyelination in a novel cohort of 50 axons. The reduced models were effective predictors; the demyelination regression (which retained less than half of the original parameters) captured 61% of the CV change observed in simulations (Figure 5C). The remyelination regression captured 87% of the observed CV recovery (Figure 5D). Comparing these parameters with empirical data, when available, may help estimate the severity of CV delays and AP failures in the cortex of aging rhesus monkeys.

Figure 5.

AP failures impair performance in a neural network model of working memory

Equipped with the quantification of impaired AP transmission due to myelin alterations in the single neuron model, we next elucidated how these impairments affected neural circuit function. We focused on spatial working memory because the underlying neural network mechanisms have been studied and modeled in detail (e.g. Compte et al., 2000; Hansel and Mato, 2013; Ibañez et al., 2020). We built on a previous spiking neural network model that accounts for many experimental findings (Hansel and Mato, 2013). It consists of 16,000 excitatory and 4,000 inhibitory integrate-and-fire neurons (see Methods). Neurons are coupled through excitatory synapses with AMPA and NMDA receptors, and inhibitory synapses with GABA receptors. Recurrent excitatory synapses are facilitating, which promotes robust and reliable persistent activity despite spatial heterogeneities in the connectivity or in the intrinsic properties of the neurons.

We first simulated the classical oculomotor delayed response task (Figure 6A) in a cohort of 10 young (control) networks with different random connectivities and intact AP transmission (see Methods). For appropriate levels of excitation and inhibition, a localized activity bump forms during the cue period, and this bump persists through the delay period. The center of the bump encodes the remembered spatial location (Figure 6B,i; Compte et al 2000; Hansel and Mato, 2013). Successful trials require a sufficiently strong activity bump throughout the delay period, quantified by the memory strength (Figure 6C, blue line). If the memory strength decreases over time (e.g. caused by the demyelination/remyelination conditions discussed below) the memory duration – the period during which the network can retain the stimulus – becomes limited (Figure 6C). Moreover, due to random fluctuations, the activity bump diffuses along the network during the delay period. This leads to trial-to-trial variability in the cue position read out from the network activity, modeling the variability of recalled spatial locations observed empirically (Figure 6B, right panels; Wimmer et al., 2014). This memory diffusion increases with the delay duration, consistent with decreasing working memory precision observed experimentally (Figure 6D; Funahashi et al., 1989). The bump movement during the delay also has a directed component, i.e., a systematic bias clockwise or counterclockwise away from the cue location, caused by heterogeneities in the network connectivity (Figure 6B, right panels and Figure 6E; Hansel and Mato, 2013). This memory drift is a possible explanation for delay-dependent systematic biases in working memory observed experimentally (see Discussion). However, because of their established relationship with working memory performance, in the following we focus on memory duration, corresponding to complete forgetting, and memory diffusion, corresponding to working memory precision.

Figure 6.

To investigate how myelin alterations affect working memory maintenance, we explored in the network model the same demyelination and remyelination conditions as we did in the single neuron model. Because our network model consists of point neurons (i.e., without detailed axons), we incorporated CV slowing as an effective increase in synaptic transmission delays (see Methods). We found that increased propagation delays had no effect on the network performance (Supplementary Figure 6). This was as expected because, by design, the network operates in an asynchronous state with irregular neural activity in which the timing of individual spikes does not affect network function. AP failures, on the other hand, did have a large impact. To simulate AP failures we created a probabilistic model of spike transmission from the excitatory presynaptic neurons to both the excitatory and inhibitory postsynaptic neurons (see Methods). For AP failure probabilities matched to the distribution of AP failure probabilities across the cohort of 50 single model neurons, we observed a decay of the activity bump (i.e., reduced memory strength over time) in the network model. This can ultimately lead to extinction of the activity bump during the delay period, representing complete forgetting of the remembered stimulus (reduced memory duration; Figure 6B, ii and iii; Figure 6C and Supplementary Figure 7A). For a mild reduction of the bump amplitude, memory duration was not affected (Figure 6C, purple line) but memory diffusion increased compared to the control network (Figure 6B, iii and Figure 6D). This increase in memory diffusion is consistent with mathematical analysis of firing rate models showing that the bump diffusion depends inversely on the squared bump amplitude (Kilpatrick and Ermentrout, 2013; Esnaola-Acebes et al., 2022).

Impact of demyelination and remyelination on working memory

We then systematically characterized changes in memory duration and memory diffusion by comparing working memory performance in the cohort of 10 control networks with the performance of those networks perturbed corresponding to varying degrees of demyelination (removing a fraction of myelin lamellae from a fraction of myelinated segments) and remyelination (replacing myelinated segments with two shorter and thinner myelin sheaths with a node in between).

Demyelination impairs working memory performance compared to control networks

We found that the memory duration was not affected when removing up to 55% of myelin lamellae per myelinated segment, regardless of the percentage of axonal segments that were altered along an axon (Figure 7A; left panel). However, when between 55% and 75% of the myelin lamellae were removed, the memory duration began to decrease, depending on the percentage of axonal segments that were demyelinated. In this case, increasing the percentage of demyelinated segments from 10% to 50% led to a progressive impairment, whereas an increase to 75% of the segments did not impair memory duration further. Finally, in cases where 100% of the myelin lamellae were removed, the memory duration dropped to ≤ 1 s, regardless of the percentage of segments that were demyelinated. In a similar trend, memory diffusion increased, that is, working memory became less precise starting when removing between 25% and 50% of the myelin lamellae (Figure 7B; left panel). Again, we found a progressive impairment, depending on the percentage of myelin lamellae removed and the percentage of myelinated segments affected, with a ceiling effect when more than 50% of segments were demyelinated.

Figure 7.

Complete remyelination recovers network function

We observed that remyelination of all previously demyelinated segments (100%), independently of the degree of demyelination, recovered the memory duration to the control-networks-like performance (Figure 7A; middle and right panels; black lines). However, working memory precision is not fully recovered in all these cases, indicated by an increase in the diffusion constant (Figure 7B). The performance of control networks was completely recovered only when 75% of the myelin lamellae were added back to the remyelinated segments. Thus, despite the new shorter and thinner myelin sheaths compared to the original intact ones, complete remyelination is able to recover control, unperturbed network function.

Incomplete remyelination leads to partial recovery

We studied the effect of incomplete remyelination (remyelination of 25% - 75% of previously demyelinated segments) after both complete and partial demyelination (Figure 7A-B, middle and right panels, respectively). When we remyelinated between 25% - 75% of the previously completely demyelinated segments, we did not observe a significant recovery of the memory duration (Figure 7A; middle panel; memory duration ≤ 1 s in all cases). Memory diffusion was partly recovered compared to the complete demyelination case (compare left and middle panels in Figures 7B) when 50% - 75% of the demyelinated segments were remyelinated, but it remained far from the performance of the control network. In sum, incomplete remyelination was unable to restore network function when bare axon (completely demyelinated) segments were present. However, when we remyelinated between 25% - 75% of the previously only partially demyelinated segments, both memory duration and memory diffusion were restored closer to the values of the control networks (Figure 7A-B; right panel).

Fewer normal myelin sheaths lead to decreased network performance

Whereas up to now we have studied network models with AP transmission probabilities corresponding to a single degree of myelin alterations, we also sought to reveal the effect on working memory performance of more biologically realistic network models with AP transmission probabilities matched to both axons with intact and with altered myelin sheaths, as likely occurs in the aging brain (Figure 1). Thus, we ran network model simulations combining AP failure probabilities corresponding to groups of neurons containing intact axons and axons presenting different degrees of demyelination (Figure 8A; Methods). This allowed us to predict working memory deficits for a degree of demyelination that is within the empirically observed range of 90% to 100% normal myelin (Peters and Sethares, 2002). We observed that the performance was impaired – memory duration significantly decreased (Figure 8B), and memory diffusion significantly increased (Figure 8C) – when the percentage of normal sheaths decreased. These results are consistent with an experimentally observed increased cognitive impairment in various learning and working memory tasks (including the delayed recognition span task, a spatial working memory task) with an age-related decrease of the percentage of normal sheaths in dlPFC of rhesus monkeys (Peters and Sethares, 2002). Importantly, our results indicate that myelin alterations alone can account for significant working memory impairment, pointing to demyelination as a key factor in age-related working memory decline.

Figure 8.

Shorter and thinner myelinated segments impair working memory

Finally, to make our model predictions regarding the effects of remyelination on working memory comparable to empirical studies of axon remyelination, we combined model neurons with intact axons and with axons containing different proportions of remyelination, and we quantified the percentage of new myelin sheaths in different groups of 50 axons (Figure 9A; Methods). We chose groups with proportions of new sheaths that matched empirical findings: percentages of paranodal profiles between 5% and 17% observed in electron microscopy images (Peters and Sethares, 2003). Increased paranodal profiles are indicative of more and shorter myelin sheaths and, thus remyelination, in the brain. We observed that performance was significantly impaired – memory duration decreased (Figure 9B), and diffusion constant increased (Figure 9C) – when the percentage of new myelin sheaths was increased. This is consistent with empirical findings showing that cognitive impairment for a cohort of 18 young and aged rhesus monkeys increased with an increase of the percentage of paranodal profiles in the dlPFC (Peters and Sethares, 2003). Unexpectedly, our model indicates that compared to the performance of networks composed of neurons possessing axons with intact myelin sheaths, both demyelination and remyelination leads to an impaired performance.

Figure 9.

Discussion

While technical limitations in EM studies make it difficult to assess the frequency of completely demyelinated segments during aging (Peters and Sethares, 2003), our single neuron model assumed that age-related myelin dystrophies – which are well-documented – alter the insulative properties of lamellae analogously to demyelination. As expected, CV in the single neuron model cohort slowed down more as demyelination increased. Once demyelinated segments lost at least half of their lamellae, several models exhibited high rates of AP failures. This finding is consistent with levels of demyelination causing CV delay and AP failure in past modeling studies (Stephanova et al., 2005; Stephanova and Daskalova, 2008; Naud and Longtin, 2019). Past modeling studies found CV slowing of 10-25% relative to unperturbed conditions after remyelinating all segments (Lasiene et al., 2008; Powers et al., 2012; Scurfield and Latimer, 2018), which is also in line with our results. An important novelty of our study is the systematic exploration of both demyelination and subsequent remyelination. We found that remyelination restored AP propagation to varying degrees, depending on severity of initial demyelination. A few model neurons had somewhat unexpected results, such as a CV increase relative to the unperturbed case when just 25% segments were completely remyelinated with 75% of lamellae (4% of the cohort). This behavior likely arises from a complex interplay between high current density imparted by ion channels in the nodes and altered electrotonic lengths in the segments where myelin changes occur (e.g., Waxman, 1978; Naud and Longtin, 2019), providing an opportunity for future study. Susceptibility to impairment or recovery from myelin dystrophy depends greatly on axonal morphology and biophysics, which we examined across our model cohort. Many parameters emerged as critical contributors to demyelination and remyelination effects, including axon diameter, node length, length of myelinated segments, and nodal ion channel densities. Parameters such as the axoplasm resistance, myelin resistance and capacitance, axon capacitance and tight junction resistance captured nonlinear interactions between the 5 morphologic parameters used to construct the cohort. Past modeling studies also found that these parameters affect axonal CV to different degrees (e.g., Goldman and Albus, 1968; Moore et al., 1978; Babbs and Shi, 2013; Young et al., 2013; Schmidt and Knösche, 2019). Better empirical measurements of these parameters in monkey dlPFC, for example from 3-dimensional EM studies, would help predict the extent to which myelin dystrophy and remyelination with aging affect AP propagation.

With the spiking neural network model, we found that increasing probabilities of AP failure, corresponding to higher degrees of demyelination, gradually impaired the time during which stimulus information can be held in working memory, as well as working memory precision. Sufficient remyelination restored performance to the control network level. In contrast to the strong impact of AP failures on network function, introducing propagation delays to mimic AP slowing did not have an effect (Supplementary Figure 6). This is because the network operates in an asynchronous state in which the dynamics are primarily governed by the statistics of neuronal activity (e.g. firing rates), rather than precise spike timings (Hansel and Mato, 2013; van Vreeswijk and Sompolinsky, 1996). While highly irregular persistent activity is indeed observed in PFC during working memory tasks (Compte et al., 2003), at the mesoscopic level oscillations of the local field potentials and synchronization also play an important role (Liebe et al., 2012; Buschman et al., 2012; Gregoriou et al., 2009; Salazar et al., 2012). AP slowing may alter these neural oscillations and synchrony which could lead to further working memory impairment not captured by our network model. In addition to age-related changes in memory duration and precision, our network model predicts an age-related increase in systematic errors due to an increased drift of the activity bump (Supplementary Figure 8). Delay-dependent systematic working memory errors have been observed in behavioral experiments (Bae, 2021; Panichello et al., 2019; Stein et al., 2021) and it would be interesting to test whether those biases also change with aging. While the prefrontal cortex is certainly a key player in working memory, it is not the only brain area involved in this function; working memory is likely sustained by the interaction of several fronto-parietal brain areas (Leavitt et al., 2017; Christophel et al., 2017). A natural extension of this work, that we will address in future studies, is the development of large-scale models of working memory (Mejías and Wang, 2022), incorporating myelin alterations in the neuronal axons forming local circuits in the cortical areas involved in and projecting to white matter tracts of inter-area connections.

For biologically realistic combinations of neurons with intact and demyelinated axons, our network model predicts that myelin dystrophies alone would lead to spatial working memory impairment with aging (Figure 8). This result supports the observation that cognitive impairment increases as the percentage of normal myelin sheaths decreases (Peters and Sethares, 2002). In addition, combining neurons containing either intact or remyelinated axons, we found that insufficient remyelination also leads to spatial working memory impairment (Figure 9). This result may explain two seemingly contradictory empirical findings. One study showed a positive correlation between cognitive impairment and a higher percentage of new, shorter and thinner, myelin sheaths (Peters and Sethares, 2003). A second study found that, with aging, oligodendroglia have a decreased capacity for effective maturation, remaining in a progenitor state without being able to produce new myelin. This low capacity for remyelination was correlated with spatial working memory impairment in the monkeys (Dimovasili et al., 2022). Our single neuron model may unify these results in the following way. First, Scurfield and Latimer (2018) found that insufficiently remyelinated axons, with large numbers of transitions between long (unperturbed) and new shorter segments, have slower CV compared to the unperturbed axons; in general we observed this too. Second, our models showed that axons with demyelinated segments, or otherwise poorly insulated myelin sheaths, experience delayed and even failed AP propagation, and in turn, lead to working memory impairment. Therefore, it is reasonable to assume that ineffective remyelination may lead to working memory impairment. In fact, complete remyelination of all the previously demyelinated segments with sufficient myelin, with fewer transitions between long and short segments, recovered working memory function. Our findings also suggest that differences in the degree of demyelination or remyelination (Figures 8 and 9) may account for the cognitive variability observed across individuals with aging, which encompass individuals with both good (successful agers) and impaired (unsuccessful agers) cognitive function (Lacreuse et al., 2005; Moore et al., 2006; Moore et al., 2017; Moss et al., 2007).

Conclusions

The multiscale modeling approach introduced here opens the door to studying the combined effect on working memory of several changes observed with aging in the monkey dlPFC (review: Luebke et al., 2010). For example, in our previous work (Ibañez et al., 2020) we modeled increased AP firing rates observed in vitro together with the loss of both excitatory and inhibitory synapses, and quantified how these alterations affected working memory performance. We found that in “aged” DRT networks that were able to maintain the memory of the stimulus, as in Wang et al. (2011), higher firing rates in individual pyramidal neurons and the loss of excitatory and inhibitory synapses partially compensated each other in order to maintain the memory. However, complete compensation only occurred when the overall excitatory drive of the network further decreased. It is conceivable that such a decrease in excitation could arise from AP failures induced by myelin dystrophy, as shown here. It has also been hypothesized that the age-related hyperexcitability of pyramidal neurons observed in vitro could work as a homeostatic regulatory mechanism to compensate for AP delays and disruptions due to dystrophic myelin (Luebke et al., 2010). The multicompartment model used here is particularly well-suited to simultaneous modeling of age-related changes in the soma, dendrites, synapses, and axon of individual pyramidal neurons. This can include changes to the nodes of Ranvier (reviewed in Arancibia-Carcamo and Attwell, 2014) and to the metabolic network associated with demyelinated neurons including higher oxygen consumption rates (Gerevich et al., 2023), which were not modeled here. To accurately quantify the effect of these three changes together on network performance, a detailed quantification of age-related myelin dystrophies and remyelination along the axons of pyramidal neurons, for example in dlPFC, would be needed.

Myelin dystrophy happens in many neurological conditions, including multiple sclerosis, schizophrenia, bipolar disorder, autism spectrum disorders, and after traumatic brain injury (Valdés-Tovar et al., 2022; Takahashi et al., 2011; Gouvêa-Junqueira et al., 2020; Franklin and ffrench-Constant, 2008; Armstrong et al., 2016; Galvez-Contreras et al., 2020; Simkins et al., 2021), in which cognitive functions are severely affected. Working memory is often one of the most vulnerable cognitive functions in these conditions, in particular in schizophrenia or autism spectrum disorders (Gold and Luck, 2023; Hahn et al., 2018; Wang et al., 2017). Thus, our findings have broader impact beyond aging and have implications for many demyelinating conditions. Overall, our study shows that myelin changes alone can account for working memory impairment, pointing to myelin degradation as a key factor in working memory decline with aging and other neurological conditions.

Methods

Single Neuron Model

To simulate age-related myelin dystrophies in individual neurons, we used a biophysically detailed multicompartment model of a rhesus monkey dlPFC Layer 3 pyramidal neuron (Rumbell et al., 2016), executed in the NEURON simulation environment (Carnevale and Hines, 2006). The soma, apical and basal dendritic arbors were constructed schematically, scaled to the overall surface area of a three-dimensional morphologic reconstruction from empirical data (68 compartments total; details in Rumbell et al. (2016)). Ion channel conductances and kinetics were fit to electrophysiological data obtained in vitro from the same neuron. In addition to passive membrane dynamics, the model included two sodium channels (fast inactivating, NaF; and non-inactivating persistent, NaP), three potassium channels (delayed rectifier, KDR; muscarinic receptor-suppressed, KM; transient inactivating A-type, KA), a high-threshold non-inactivating calcium channel (CaL), a calcium-dependent slow potassium channel (KAHP), and the hyperpolarization-activated anomalous rectifier channel (AR). Because our prime focus was quantifying alterations in action potential (AP) propagation along the axon under dystrophic myelin conditions, we augmented the axon hillock and initial segment of the Rumbell et al. (2016) model (5 compartments each) by attaching nodes and myelinated segments from a detailed axon model (Gow and Devaux, 2008; Scurfield and Latimer, 2018). To the initial segment we attached 101 nodes (13 compartments each) alternating with 100 myelinated segments. Each myelinated segment was bound by a group of four paranodes (5 compartments each) on either side, with an internode flanked by two juxtaparanodes (9 and 5 compartments each, respectively) sandwiched between the paranode groups (Figure 2A). On both ends of a myelinated segment the extreme outward paranode interfaced with the adjoining node. Segments between successive nodes were endowed with myelin lamellae (wraps), and with tight junctions between the innermost lamella and the axolemma necessary for improving insulation and accurate modeling of AP propagation in nerve fibers with diameters less than 0.9 μm as often observed in PFC. Axon compartments included passive membrane dynamics as well as the NaF and KDR channels in the nodes. Simulations used a fixed 0.025 ms time step.

As done during in vitro electrophysiological experiments (Chang et al., 2005; Ibanez et al., 2020) and past modeling studies (Coskren et al., 2015; Rumbell et al., 2016), we first applied a holding current to stabilize the somatic membrane potential at -70 mV, then injected a +380 pA current step into the somatic compartment for 2 seconds. We recorded the somatic firing rate, as well as the conduction velocity (CV) of APs propagating from the first node after the initial segment (proximal to the soma) and the penultimate node (at the distal end). The CV was relatively insensitive to variations in the magnitude of suprathreshold somatic current steps, and whether the current was constant or included Gaussian noise. Therefore, we quantified CV changes and AP failures from responses to constant +380 pA current steps.

Building a cohort of control model neurons

Recognizing that the effects of myelin dystrophy may depend on the physical dimensions of an axon, we constructed a cohort of neuron models that were consistent with empirical measurements of axonal morphology in young adult rhesus monkeys. To form this ‘control’ model cohort, we identified five parameters of axonal morphology which have been measured empirically in rhesus monkey cortex: axon diameter (Peters et al., 2001; Bowley et al., 2010; D.L. Rosene, unpublished observations), node length (Peters and Sethares, 2003), length of myelinated segments (Waxman, 1980), number of myelin lamellae (Peters and Sethares, 2003), and thickness of each lamella (Moore et al., 1978; Peters and Sethares, 2003). We also defined three scale factors defining the ratio of leak, NaF, and KDR conductances in our nodes relative to those of the Scurfield and Latimer (2018) model, assuming an upper bound of 1 for each. This gave a total of eight parameters with associated feasible parameter ranges that we varied (Table 1); other parameters including tight junction resistivity, myelin resistivity, and myelin capacitance were held fixed at values from Scurfield and Latimer (2018). We then used the space-filling Latin Hypercube Sampling (LHS) design (Morris and Mitchell, 1995; Rumbell et al., 2016; Ibañez et al., 2020) to identify a set of 1600 points in parameter space that maximized the minimum distance between all pairs of points. We simulated each model specified by these points, and had several selection criteria for identifying biologically plausible models. First we constrained somatic firing rates within ranges observed empirically in Layer 3 pyramidal neurons of rhesus monkey dlPFC: firing 13-16 Hz in response to the +380 pA current step and silent when no current was injected (Chang et al., 2005; Ibañez et al., 2020). We also required CVs of 0.3 –0.8 m/s (rhesus monkey corpus callosum, D.L. Rosene, unpublished observations), and ensured that simulated APs were suprathreshold in the nodes and subthreshold in the juxtaparanode and internode regions, indicating saltatory conduction. Of the 1600 simulated models, 138 met these criteria; for the present study, we randomly selected 50 models to comprise the control model cohort. The g-ratio (ratio of axon to fiber diameter) among models in the cohort was 0.71 ± 0.02, with total axon lengths of 1.2 ± 0.1 cm.

Table 1

Simulating myelin alterations

We assumed that the effect of all myelin dystrophies observed empirically (e.g. Figure 1) could be modeled by removing lamellae from myelinated segments (demyelination, which reduces electrical insulation), and that remyelination could be modeled by replacing myelinated segments with two shorter and thinner segments with a node in between. To simulate demyelination, we varied two independent factors: the percentage of myelinated segments selected for demyelination along an axon (demyelination percentage), and the percentage of myelin lamellae removed from those segments (lamellae removal percentage). For each demyelination percentage (10, 25, 50, and 75%; Figure 3), we generated 30 randomized lists of segments to demyelinate. Then for each of the 30 trials, each lamellae removal percentage was applied (25, 50, 55, 60, 65, 70, 75, and 100%) to the chosen segments, for all 50 models in the cohort. We then simulated the +380 pA current step, calculating the CV and the number of APs that propagated to the distal end of the axon. For each perturbation, we defined the CV change as the percentage change in CV induced relative to the CV of the corresponding unperturbed model. We also computed, for each perturbation, the percentage of AP failures at the distal axon end, relative to the number of APs observed at the first node. Trials with AP failure rates of 75% or more were labeled as ‘high failure trials’.

To simulate remyelination three factors were varied: the percentage of myelinated segments initially demyelinated (demyelination percentage: 25, 50, 75%); the percentage of those affected segments which were then remyelinated (remyelination percentage: 25, 50, 75, 100%); and the percentage of lamellae restored with remyelination (lamellae restoration percentage: 10, 25, 50, 75%). We performed these remyelination simulations under two demyelination conditions: where the initially demyelinated segments had lost all their lamellae (‘complete demyelination’), or had lost half their lamellae (‘partial demyelination’). Sample remyelination perturbations shown in Figure 4. Remyelination was performed by replacing an affected (previously demyelinated) segment with two shorter segments, each including paranodes, juxtaparanodes, and an internode, and a new node between them. The number of lamellae on the shorter segments was determined by the lamellae restoration percentage. As before, for each demyelination percentage (25, 50, and 75%), we generated 30 randomized lists of segments to demyelinate. Segments to remyelinate were decided based on remyelination percentage. For example, in the case of 50% remyelination every alternate (previously) demyelinated segment was remyelinated. Then for each trial, we applied each combination of remyelination percentage and lamellae restoration percentage to the chosen segments, for all models in the cohort; and simulated the current clamp protocol. For each remyelination condition, in addition to computing the percentage of high failure trials, we defined CV recovery as the percentage improvement in CV relative to the CV change for the completely demyelinated case. For example, if a control model had a CV of 1 m/s, and a CV of 0.6 m/s after complete demyelination, then the CV change was -0.4/1.0 = -40%. If a subsequent remyelination condition led to a CV of 0.8 m/s (an increase of 0.2 m/s over the demyelinated case), the CV recovery was 0.2/0.4 = 50%.

Statistical assessment of model parameter importance

To identify which parameters of the multicompartment model had the greatest influence on axonal responses to demyelination and remyelination perturbations, we used least absolute shrinkage and selection operator (lasso) regression (Tibshirani, 1996; James et al., 2021), implemented in MATLAB R2022a (Mathworks, Natick, MA). Lasso fits a regression model to response variable, Y = (y₁, y₂, . . ., y_n)^T, given predictor variables X (n observations of p predictors) by selecting coefficients that minimize the quantity

for a given value of λ (James et al., 2021). The former quantity is the residual sum of squares of the model. The latter quantity is the absolute sum of coefficients; including it in the minimization shrinks some coefficients to zero when λ is sufficiently large, enabling lasso to perform feature selection. The tuning parameter λ was selected to minimize the 10-fold cross validation error; then for the selected value of λ, the regression was repeated using all available data. A total of 12 parameters served as predictor variables for the lasso regression. We started with 6 of the 8 parameters used to construct the cohort (axon diameter, node length, myelin length, and scale factors for leak, NaF, and KDR conductances), and combined two others parameters (number of myelin lamellae and lamella thickness) into one quantity: myelin thickness (their product). Past studies (Goldman and Albus, 1968; Koles and Rasminsky, 1972; Moore et al. 1978; Gow and Devaux, 2008) identified several other parameters that affect axonal propagation, including the g-ratio, axoplasmic resistance, axon capacitance, myelin resistance, myelin capacitance, and tight junction resistance. We added five of these parameters to the seven cohort parameters as predictor variables, omitting only g-ratio due to its high correlation with axon diameter and myelin thickness. We created response variables summarizing the effects of demyelination and remyelination on CV in each of the 50 members of the model cohort. AP failures were recorded as CV change of -100%. We did not include high failure trials as a response variable, since CV changes precede AP failures. For demyelination, we averaged the CV change for all randomized trials in which 100% of lamellae were removed from 25%, 50% and 75% of segments. For remyelination, we averaged the CV recovery in all randomized trials in which 25, 50, or 75% of segments were completely demyelinated, and then all affected segments were remyelinated as two shorter segments with 75% of lamellae added back. To facilitate comparison, we z-scored all predictor and response variables before performing lasso separately on each response. We tested the predictive ability of the lasso models by randomly selecting another 50 of the 138 models from the original hypercube that met the inclusion criteria, then simulating the specific myelin alteration protocols that comprised the demyelination and remyelination response variables. We computed Pearson’s correlation coefficient between the z-scored predicted vs. observed responses.

Spiking neural network model

We adapted a neural network model (Hansel and Mato, 2013) to simulate the neural circuit in the dlPFC underlying spatial working memory during the oculomotor delayed response task (Figure 6A). The network model is composed of N = 20,000 leaky integrate-and-fire neurons, N_E = 16,000 excitatory (E) neurons (80%) and N_I = 4,000 inhibitory (I) neurons (20%) with sparse probabilistic connections among all neuronal populations. A full description of this probabilistic version of the ring model can be found in Hansel and Mato, 2013. Here we describe the essential features of the network and summarize all the parameter values in Supplementary Table 1. The range of the interactions is represented by the parameters σ_EE, σ_EI, σ_IE, and σ_II. Each neuron receives K = 500 total inputs, where K_E = 0.8 K are excitatory inputs and K_I = 0.2 K are inhibitory inputs. The subthreshold membrane potential of each excitatory and inhibitory neuron (i) in the network is described by

where τ_E and τ_Iare the membrane time constants for the E and I neurons, respectively. is the total recurrent synaptic current that each neuron receives from all the other neurons in the network connected to it. is a constant background input, representing internal brain currents that come from outside the network. represents the transient sensory input to each E neuron, associated to a given direction of the stimulus, and only active during the cue period of the task. An action potential is fired each time that the membrane potential of a neuron reaches the threshold value, V_T. The voltage of the membrane is reset to the baseline value, V_R, immediately after.

The total recurrent synaptic current for each neuron, is given by

where subindices a and n represent the excitatory AMPA and NMDA glutamatergic receptors, and g an inhibitory GABA receptor. Each synaptic current is given as in (Hansel and Mato, 2013), with synaptic decay time constants for each receptor, τ_a, τ_n, and τ_g, respectively. Short-term plasticity was incorporated to the excitatory-to-excitatory synaptic connections through the variables u and x, given by (Markram et al., 1998)

The sensory input current to excitatory neuron i during the cue period, and for a specific location of the stimulus, is given by.

where θ_i is the preferred direction of the neuron i, given by its position on the ring, θ_cue is the cue direction, and I_max,E and ε_E are the amplitude and the width of the sensory input current, respectively.

Simulating a cohort of young (control) networks

Using the Brian2 simulator (based on Python), we simulated 10 different networks, all with the same parameters, but each with a different connectivity profile. For each network, we ran 280 trials, with different initial conditions and cue positions. All networks performed the DRT – with 2 seconds fixation period, 1 second cue period and 4 seconds delay period – and maintained the memory of the stimulus during the whole delay period of the task (Figure 6B,i and Supplementary Figure 5A, left panel).

Modeling the effects of demyelination and remyelination in the network model

Modeling action potential propagation failures in the network

To model the action potential failures at the distal end of the axons quantified with the single neuron model under the different demyelination and remyelination conditions, we perturbed the 10 control networks by designing a probabilistic model of spike transmission from the excitatory presynaptic neurons to both the excitatory and inhibitory postsynaptic neurons. From the single neuron model, for each demyelination/remyelination condition, we quantified the probability of AP failure for each of the neurons in the control cohort, as well as the percentage of those neurons that shared the same probabilities of failure. That is, the percentage of neurons that had probability of failure = 0, probability of failure = 1 or any other probability. Then, we computed the probability of transmission, P_transmission = 1 − P_failure, and we specified P_transmission for the corresponding percentages of excitatory neurons in the networks. Thus, in the network model we took into account the heterogeneity observed in the single neuron model for the different demyelination/remyelination conditions.

Modeling conduction velocity slowing in the network

To explore the effect of CV slowing along the axons of model neurons, we simulated 20 control networks and 20 networks for the case of demyelination of 50% of the segments along the axons of model neurons by removing 60% of the myelin lamellae (we ran 280 trials for each network). Then, we added random delays – with a minimum value of 0 ms and a maximum value of 100 ms in the control networks, and of 40 ms in the perturbed networks – in both the AMPA and NMDA excitatory connections to both E and I neurons (Supplementary Figure 6).

Quantification of normal myelin sheaths in groups of model neurons containing both intact and demyelinated axons

We created different groups of 50 total model neurons containing, among those 50, different random amounts of neurons with intact axons and with demyelinated axons, with either 10, 25, 50 or 75% of segments demyelinated. In each group, the distribution of the demyelinated segments along the altered axons was also randomly chosen among the 30 possible distributions simulated with the single neuron model. We sorted the axons in each group by locating their origin aligned in the same position, divided the axons longitudinally in sections of 0.5 μm, and calculated the percentage of normal, unperturbed myelin sheaths in each section (Figure 8A). Then, we averaged across all sections and we picked 60 groups that had over 80% of normal myelin sheaths. According to the proportion of intact and demyelinated axons (with either 10, 25, 50 or 75% of segments affected) in each group, we adjusted the distribution of AP transmission probability (P_transmission; P_transmission = 1 for intact axons) in 1 of the 10 control networks. To do this we considered that the perturbed axons in 40 of the groups were completely demyelinated (all lamellae removed) and that in the remaining 20 groups, they had 75% of lamellae removed.

Quantification of new myelin sheaths in groups of model neurons containing both intact and remyelinated axons

We again created different groups of 50 total model neurons containing different random amounts of neurons with intact and remyelinated axons. The remyelinated axons had either 25, 50 or 75% of segments remyelinated (with two shorter and thinner myelin sheaths), following previous demyelination of either 25, 50 or 75% of the segments. In each group, the distribution of the remyelinated segments along the altered axons was randomly chosen among the 30 possible distributions. As before, we sorted the axons in each group by locating their origin aligned in the same position, we divided the axons longitudinally in sections of 0.5 μm, and calculated the percentage of remyelinated segments in each section (Figure 9A). Then, we averaged across all sections, we multiplied by two to calculate the number of new myelin sheaths, and we picked 22 groups that had below 45% of new myelin sheaths. According to the proportion of intact and remyelinated axons (with different percentages of remyelinated segments) in each group, we adjusted the distribution of P_transmission in the same control network as before for the case of remyelination following partial demyelination by adding 25% of the lamellae back.

Measures of working memory performance

The remembered location was obtained using a population vector decoder and network performance was quantified by the following four measures that describe the quality of the activity bump maintaining a memory of the stimulus during the delay period: the memory strength, the memory duration, the drift rate and the diffusion constant.

Memory strength

The memory strength was defined as the modulus M(t) of the population vector, Z(t), which characterizes the spatial modulation of the excitatory neuronal activity at time t, given by

Decoded cue location

Ψ(t) is the argument of Z(t) and indicates the remembered stimulus location. That is, the center of the activity bump (Figure 6B).

Memory duration

The memory duration was defined as the time from the delay onset until the time where the memory strength decayed below a fixed threshold that we set at 0.4. That is, the memory duration is the interval during which the bump of neural activity is reliably maintained (Figure 6C).

Memory drift

The drift was defined as the bias of the estimator, which is given by

That is, the difference of the trial-average of the estimates Φ(t) and the true value θ_cue (Esnaola-Acebes et al., 2022; Figure 6E). The drift rate is the slope of a linear fit of the drift during the memory duration time (Supplementary Figure 8).

Memory diffusion

The diffusion was defined as the variance of the estimator, given by

the variation of the estimates about their mean value (Esnaola-Acebes et al., 2022; Figure 6D). The diffusion rate is the slope of a linear fit of the diffusion during the memory duration time (Figure 7B).

Acknowledgements

We acknowledge the use of Fenix Infrastructure resources, which are partially funded from the European Union’s Horizon 2020 research and innovation programme through the ICEI project under the grant agreement No. 800858, and the research cluster at Franklin and Marshall College, funded through NSF grant 1925192. We thank David Latimer and Albert Compte for sharing computer code, as well as Jason Brooks and Tony Weaver for technical assistance. We thank the CERCA Programme/Generalitat de Catalunya for institutional support.

Funding

This work was supported by NIH/NIA grant R01 AG059028, NIH 1R01AG071230-01, and grant PCI2020-112035 from MCIN/AEI/10.13039/501100011033 and the European Union “NextGenerationEU”/PRTR. This work was supported by the Spanish State Research Agency, through the Severo Ochoa and María de Maeztu Program for Centers and Units of Excellence in R&D (CEX2020-001084-M).

Software and data availability

The single neuron and network models for this study will be made available on ModelDB upon publication.