All simulations were performed using NEURON, which is embedded in Python 2.10 (Hines et al., 2009). The backward Euler method was used for numerical integration with a time step of 0.025 ms. Two full-morphology compartmental models of rat CA1 pyramidal cells (Kim et al., 2015; Magó et al., 2020) were used for the simulations of the LTP experiments. The simulation files of the Kim et al. (2015) CA1 model were downloaded from ModelDB (https://modeldb.science/), accession number 184054. We added 150 excitatory synapses, randomly distributed on the apical tuft dendrites of the model (Fig. 2a). While spines were not explicitly modelled, the associated surface area was accounted for by adjusting the specific membrane resistivity (Rm) and specific membrane capacitance (Cm) of compartments located more than 100 μm from the soma, which were multiplied by a factor of two. Each synapse was composed of an AMPA and an NMDA conductance, simulated by the sum of two exponential functions with rise time and decay time constants of 0.2 and 2 ms for AMPA (Katz et al., 2009) and 1 and 50 ms for NMDA (Spruston et al., 1995). Initial peak conductances were randomly selected from a lognormal distribution (mean 0.18, sigma 0.35 nS) for both AMPA and NMDA synapses (Rößler et al., 2023). The voltage-dependent magnesium block of the NMDAR was simulated using the equation: \(_=_^\times exp\left(-0.062\times \left(V-10\right)\right)\right]}^\), where \(_^=1\) mM is the Mg2+ concentration in the bath and \(V\) is the local dendritic voltage (Kim et al., 2015). In addition, we decided to modify the model by turning off the slow inactivation of Nav channels, which resulted in a better fit to the experimental data. This modification is supported by a review of the literature and published models of CA1 cells (Bloss et al., 2018; Jarsky et al., 2005), as the presence or absence of slow inactivation in sodium channels varies across different neuronal models and experimental observations.

Simulations of suprathreshold LTP in apical tuft dendrites of the CA1 pyramidal cell. (a) Random placement of 150 excitatory synapses with AMPA and NMDA receptors along the apical tuft dendrites. (b) Comparison of the plasticity results obtained from the simulations of four different theta-burst stimulation (TBS) protocols with experimental results from Kim et al. (2015). The average LTP 1 min after TBS is shown as the mean of all synaptic weights ± SEM for each protocol. (c) Detailed analysis of the first burst under different stimulation protocols. Synaptic weight changes are shown as a function of time. (d) Effect of TTX on Na-dSpike generation. Sodium currents were measured from the same dendritic location in both control and TTX conditions. (e) Model predictions for LTP as a function of distance from the soma for different stimulation protocols, including those under TTX conditions

The model of Magó et al. (2020) from the ModelDB (accession number 265511) included synaptic conductances: AMPA had a rise time of 0.1 ms, a decay time of 1 ms, and a maximum conductance of 0.6 nS, while NMDA had a rise time of 2 ms, a decay time of 50 ms, and a maximum conductance of 0.8 nS. The voltage-dependent magnesium block was simulated using the equation \(_=__^/4.3\times exp\left(-0.071\times V\right)\right)}^\) where \(_^=1\) mM is the Mg2+ concentration in the bath and \(V\) is the local dendritic voltage (Magó et al., 2020). In this model, synapses were placed on high-impedance dendritic spines consisting of a spine neck (length: 1.58 μm; diameter: 0.077 μm) and a spine head (length: 0.5 μm; diameter: 0.5 μm) with a total neck resistance of ~ 500MΩ (Harnett et al., 2012). Two, three, four or eight spines were placed on distal dendritic segments of 5 selected perisomatic dendrites (x = 0.96, Fig. 3a). To account for spines, Cm was increased, and Rm was decreased by a factor of 2 in dendritic compartments beyond 100 μm from the soma (Magó et al., 2020).

Simulations of the subthreshold and suprathreshold LTP in the perisomatic dendrites of the CA1 pyramidal cell. (a) Location of high-impedance dendritic spines with activated synapses (red dots) at the distal ends of highlighted perisomatic dendrites (dark blue). (b) Summary plot of induced LTP compared with experimental data from Magó et al. (2020) (mean ± SEM). (c) Evolution of synaptic weights of individual synapses for different scenarios. (d) Voltage traces from one spine head. (e) Dendritic voltage (top panel) with the − 37 mV threshold horizontal line and the evolution of the corresponding synaptic weights during the entire stimulation protocol

To model synaptic plasticity, the AMPA conductance representing synaptic weight was modified according to the ETDP synaptic plasticity rule (Fig. 1). In the ETDP, a presynaptic event is a presynaptic spike, while a postsynaptic event is detected when the voltage at the synaptic site exceeds the threshold of -37 mV (Jedlicka et al., 2015). Nearest-neighbor pairing was used to match presynaptic and postsynaptic events, where each presynaptic event was paired with a postsynaptic event that occurred before it and one after it. Weight change is calculated using the formula: \(w\left(t+\delta t\right)=w\left(t\right)\times \left(1+_-_\right)\), where \(_\) is positive weight change (potentiation, LTP) and \(_\) is negative weight change (depression, LTD). Potentiation occurs when the presynaptic event preceded the postsynaptic event. Conversely, depression occurs when the postsynaptic event precedes the presynaptic spike. Thus, \(_\) and \(_\) are calculated according to the formula: \(_\left(t\right)=_\times exp\left(-t/_\right)\) if \(t>0\) and \(_\left(t\right)=_\times exp\left(t/_\right)\) if \(t<0\), where \(t=_-_\), \(_\) and \(_\) are potentiation and depression amplitudes, respectively, \(_\) and \(_\) are decay constants for the time windows of plasticity. The values of \(_=0.009\)\(_=0.0012\) for the TBS protocol, \(_=0.0035\), \(_=0.001\) for the LFS protocol and \(_=_=15 ms\) for both protocols were optimized by hand. The model is very sensitive to the value of the postsynaptic event threshold, while LTP and LTD amplitudes and decay constants only influence the quantitative match.

For suprathreshold LTP, the 2stim_3xTBS stimulation protocol consisted of 3 trains of 2 pulses delivered at 100 Hz with a theta (5 Hz) interburst frequency repeated 3 times at 4 s intervals. The 5stim_3xTBS stimulation protocol consisted of 3 trains of 5 pulses delivered at 100 Hz with a theta (5 Hz) interburst frequency, repeated 3 times at 4 s intervals. Each 5stim_3xTBS protocol was simulated: (1) paired with brief (2 ms) somatic current injections at 50 Hz to elicit 3 action potentials during each burst (5stim_3xTBS_IClamp), (2) with the soma voltage clamped at − 70 mV (5stim_3xTBS_VClamp), or (3) alone (5stim_3xTBS) (Kim et al., 2015). The protocol used to induce subthreshold and suprathreshold LTP in perisomatic dendrites consisted of delivering 50 quasi-synchronous stimulations of selected spines (with a stimulus interval of 0.1 ms between spines) at a frequency of 3 Hz. Before and after the stimulation protocol, a set of four spines was stimulated separately, with a 200 ms interval between spines, and the trials were repeated at 0.5 Hz (Magó et al., 2020). Each simulation involved clustered spines on a single distal dendritic segment, with multiple simulations conducted for different dendrites to include the impact of spine positioning on synaptic plasticity.