We constructed a mathematical representation of synaptic population dynamics, modeling synapses existing in four discrete states: (1) Adult synapses with GluN2A-dominant NMDA receptors, (2) Juvenile synapses with GluN2B-enriched receptors, (3) Silent synapses lacking functional AMPA receptors, and (4) Matured silent synapses with calcium-permeable AMPA receptors (Figure 1A). This discrete-state approach, while simplified compared to the continuous spectrum of synaptic properties observed biologically, allows for tractable mathematical analysis of population-level transitions.

Critically, our model incorporates two distinct biological processes:

Process 1 - Rejuvenation: Receptor composition switching in existing synapses (Adult ↔ Juvenile), representing GluN2A-to-GluN2B subunit changes without synapse creation or elimination. This process conserves the number of functional synapses: N_adult(t) + N_juvenile(t) = N₀ at all times (where N₀ = 1000 is the baseline count), providing an internal consistency check.

Process 2 - Synaptogenesis: De novo generation of silent synapses during drug exposure, consistent with experimental observations that cocaine exposure generates new silent synapses that add to (rather than replace) the baseline neural circuitry (Huang et al., 2009; Russo et al., 2010). This is a key distinction from models that assume fixed total synapse counts.

Unlike models that assume a fixed synapse count, our framework explicitly tracks a dynamic total synapse population:

This total is not conserved. Rather:

Mathematically, summing the differential Equations 3–6 yields:

This confirms that total synapse count changes due to genesis (positive, during exposure) and pruning (negative, during withdrawal), consistent with experimental observations (Huang et al., 2009; Russo et al., 2010).

The drug exposure function D(t) is defined as a binary pulse train representing intermittent cocaine exposure sessions:

where t_start = 100 time units is the first exposure onset, ΔT = 30 time units is the inter-exposure interval, τ = 5 time units is the duration of each exposure session, and the protocol consists of 5 intermittent exposures. Under our time scaling (1 time unit = 2 h; see Section 2.6), this corresponds to 5 sessions of ~10 h each, separated by ~2.5 days. This intermittent extended-access style protocol is a stylized representation of repeated cocaine exposure schedules used in self-administration and withdrawal studies, where prolonged access can produce stronger and more persistent neuroadaptations than short-access regimens (Ahmed and Koob, 1998; Zimmer et al., 2012).

The model dynamics are governed by the following differential equations (Figure 1C). A key feature is that state transitions are explicitly phase-gated to ensure biological consistency.

Process 1 - Rejuvenation dynamics (receptor composition switching):

Note that Equations 4, 5 sum to zero, ensuring conservation: ddt(Nadult+Njuvenile)=0.

Biological rationale for directionality: During drug exposure (D(t) = 1), cocaine-induced signaling drives GluN2A-to-GluN2B receptor subunit switching (Adult → Juvenile), consistent with experimental observations of rapid GluN2B upregulation following cocaine (Huang et al., 2009; Wang et al., 2021). During withdrawal [D(t) = 0], homeostatic mechanisms promote gradual recovery toward adult-like receptor compositions (Juvenile → Adult), consistent with observed partial normalization during abstinence (Dong and Nestler, 2014).

Process 2 - Silent synapse dynamics (de novo generation and phase-gated fate):

Several important features of Equations 6, 7:

The explicit phase-gating ensures that the mathematical formulation is consistent with the biological narrative:

During drug exposure (D(t) = 1):

During withdrawal (D(t) = 0):

To facilitate interpretation, we define explicit time scaling: 1 time unit = 2 h. This mapping is derived from experimental time courses reported in the literature:

Under this scaling, the simulation protocol corresponds to:

Rate constants were selected to match relative time scales suggested by experimental observations (Table 1). Below we provide the experimental rationale for each parameter:

To explore the functional consequences of population changes, we defined theoretical indices for plasticity capacity and memory strength (Figure 1D). Total plasticity capacity was calculated as a weighted sum:

where N₀ = 1, 000 is the baseline synapse count used for normalization. The weights reflect hypothetical relative plasticity capacities based on known properties of different receptor compositions (Table 1). Note that the weight w_mature = 3.0 reflects the potentiation efficacy of CP-AMPAR-containing synapses—their ability to amplify LTP expression due to high single-channel conductance—rather than metaplasticity or ease of modification.

The theoretical memory strength index evolved according to distinct mechanisms in each phase, with saturation to prevent unbounded growth:

where α = 0.5, β = 0.1, and M_max = 30 is the saturation level.

Rationale for separate equations: The distinct mechanisms during exposure versus withdrawal reflect fundamentally different biological processes:

Saturation prevents unbounded growth: The logistic factor (1 − M/M_max) ensures that memory strength approaches an asymptote rather than growing indefinitely. This is biologically realistic: craving incubation eventually plateaus in experimental studies (Pickens et al., 2011).

NMDA receptor composition was tracked to include GluN2B content across multiple synapse types:

This formulation explicitly accounts for GluN2B presence in different synapse types:

To assess model robustness and address concerns about deterministic model stability, we performed systematic parameter sensitivity analysis by varying each key parameter (k_{a → j}, k_genesis, k_maturation) by ±50% from baseline values. For each parameter combination, we quantified effects on three output metrics: final memory strength, peak juvenile synapse count, and final mature synapse count. This analysis demonstrates that qualitative predictions are preserved across a wide range of parameter values.

Initial condition robustness: Initial conditions were fixed to represent the experimentally motivated baseline state (adult synapses only, N₀ = 1, 000, with N_juvenile = N_silent = N_mature = 0). We additionally tested alternative initial states with ±10%–20% juvenile synapses at baseline and confirmed that qualitative conclusions (coordinated population transitions, biphasic memory strengthening, and incubation dynamics) remain unchanged. This robustness reflects the fact that model dynamics are driven primarily by the drug exposure protocol rather than initial state perturbations.

To address the specificity of neural rejuvenation to addictive drugs versus natural rewards, we performed comparative simulations. Natural rewards were modeled with substantially reduced rejuvenation parameters: k_genesis = 0 (no silent synapse generation) and k_{a → j} = 0.008 (10% of cocaine effect). This parameterization is based on experimental observations that:

We note that some natural reward paradigms may show incubation-like effects under certain conditions; our comparison specifically models the absence of the silent synapse generation mechanism.

Several key assumptions underlie our model framework:

All differential equations were solved numerically using Euler's method with a time step of 0.1 time units (corresponding to 12 min biological time), implemented in MATLAB R2023a. The simulation was run deterministically. Parameter sensitivity analysis involved 15 additional simulations (5 variations × 3 parameters). Initial condition robustness was verified with 4 additional simulations testing ±10% and ±20% juvenile synapses at baseline. Code is available at the repository listed in the Data Availability Statement.

The computational framework demonstrates coordinated reorganization of synaptic populations that qualitatively parallels experimental observations of neural rejuvenation (Figures 2A–F). At baseline, the entire theoretical population consisted of adult-type synapses (N₀ = 1, 000 synapses). The simulated drug exposure protocol involved five intermittent exposures between time units 100–220, followed by extended withdrawal until time unit 500.

Upon simulated drug exposure initiation, rapid conversion to juvenile-like states occurred in the model, with the juvenile population increasing substantially during exposure periods, reaching peak levels of approximately 500 synapses. This Adult → Juvenile conversion reflects cocaine-induced GluN2A-to-GluN2B receptor subunit switching. Concurrent with this conversion, de novo silent synapse generation emerged during exposure periods, reaching maximum populations of approximately 400 synapses (constrained by the carrying capacity K_max = 500).

Figures 3A–F summarizes the model's silent synapse dynamics, incubation-related measures, experimental predictions, and overall synapse population changes across phases. Importantly, the total synapse population increased dynamically from the baseline of 1,000 to approximately 1,400 synapses during the exposure phase (Figure 3F), reflecting de novo synaptogenesis. This 40% increase is consistent with experimental observations (Huang et al., 2009). During withdrawal, the total decreased to approximately 1,300 synapses due to pruning of unmaturated silent synapses, but remained elevated above baseline due to the persistent population of matured synapses.

During withdrawal, the juvenile population gradually declined (Juvenile → Adult recovery via homeostatic mechanisms) while a significant population of matured silent synapses emerged and stabilized at approximately 300 synapses by simulation end. Critically, due to the phase-gated equations, maturation and pruning occurred only during withdrawal periods, consistent with experimental observations (Lee et al., 2013).

The simulated synaptic population transformation corresponded to profound changes in theoretical NMDA receptor subunit composition (Figure 2C). Under baseline conditions, the model exhibited mature characteristics with GluN2A receptors comprising approximately 80% of the population and GluN2B receptors representing approximately 20%, consistent with adult synaptic states.

Simulated drug exposure triggered dramatic reversal of this pattern. GluN2B-containing receptors increased rapidly during exposure periods, reaching peak levels of approximately 90% (reflecting contributions from juvenile synapses, GluN2B-rich silent synapses, and early mature synapses that retain some GluN2B), while GluN2A receptors declined correspondingly. During withdrawal, gradual recovery toward adult-like ratios occurred, though GluN2B levels remained elevated at approximately 40% compared to baseline by simulation end.

The model's memory formation index demonstrated continuous enhancement throughout exposure and withdrawal phases via distinct theoretical mechanisms (Figures 2D, 3B). During simulated drug exposure, memory strength increased in discrete increments corresponding to individual exposure episodes, driven by the enhanced plasticity capacity of the rejuvenated circuit. Memory strength increased from baseline (0) to approximately 10 units by the end of the exposure phase.

During withdrawal periods, memory strength exhibited sustained growth despite absence of drug exposure, driven by the ongoing silent synapse maturation flux (Figure 3E). The memory formation rate during withdrawal depends on k_maturation·N_silent—the rate at which silent synapses are actively recruiting CP-AMPARs. As the silent synapse pool is depleted (through maturation and pruning), this flux decreases, naturally producing the characteristic “deceleration” of incubation observed experimentally (Pickens et al., 2011). The incubation index demonstrated approximately 18-fold enhancement relative to baseline, eventually approaching the saturation level M_max.

Population-level changes translated into significant alterations in the model's overall plasticity capacity index (Figure 2E). Baseline plasticity was normalized to 1.0. Simulated drug exposure rapidly increased total plasticity capacity to peak levels of approximately 2.0–2.5 during active exposure periods, reflecting the enhanced potentiation efficacy of juvenile synapses (GluN2B-mediated Ca²⁺ dynamics) and accumulating mature synapses (CP-AMPAR high conductance).

During withdrawal, plasticity capacity remained elevated at approximately 1.8–2.0 above baseline due to persistent populations of matured silent synapses. This sustained elevation illustrates how rejuvenation effects might theoretically persist long after drug exposure cessation, creating a prolonged “vulnerability window” for relapse. This is also reflected in the increased LTP capacity (Figure 3D).

Silent synapse populations underwent distinct temporal evolution that provides mechanistic insight into incubation phenomena (Grimm et al., 2001; Kalivas and O'Brien, 2008) (Figure 3A). NMDA-only silent synapses were rapidly generated de novo during simulated drug exposure, reaching peak levels of approximately 400 synapses. During withdrawal, this population declined through competitive processes of maturatio and elimination—processes that are explicitly phase-gated to occur only during withdrawal.

Matured silent synapses emerged specifically during withdrawal, representing successful CP-AMPAR recruitment. This population increased progressively, stabilizing at approximately 300 synapses by simulation end. The rate of this maturation process (the flux) drives memory strengthening during withdrawal, explaining why craving can continue to increase even as the rate of new mature synapse accumulation slows.

Systematic parameter sensitivity analysis demonstrated that the model's qualitative predictions are robust to parameter variations (Figure 4A). Varying each key parameter by ±50% produced the following findings:

Initial condition robustness: We additionally tested alternative initial states with ±10–20% juvenile synapses present at baseline (instead of 100% adult synapses). All qualitative conclusions remained unchanged, confirming that model dynamics are driven primarily by the drug exposure protocol rather than initial state perturbations.

These results demonstrate that the deterministic model produces stable, reproducible conclusions across a wide range of biologically plausible parameter values and initial conditions.

Comparative simulations between cocaine and natural reward exposure revealed striking differences in neural rejuvenation dynamics (Figures 4C, D). Cocaine exposure produced:

In contrast, natural reward exposure (modeled with k_genesis = 0) showed:

These predictions are consistent with experimental observations that natural rewards do not recruit the silent-synapse/NMDAR rejuvenation signature seen with cocaine (Huang et al., 2009).

The model generates phase-specific experimental predictions (Figure 3C) that can guide future investigations. The bar chart shows end-of-phase snapshot values for GluN2B enrichment and total silent synapse populations:

Building on Nestler's foundational neural rejuvenation hypothesis (Dong and Nestler, 2014; Nestler, 2001), our computational approach addresses a key gap in current understanding: the population-level coordination of rejuvenation processes and their collective impact on circuit function.

Our model directly implements Nestler's core mechanisms within a population dynamics framework, with critical distinctions:

A key contribution of our model is the mechanistic explanation for continued memory/craving growth during withdrawal. The memory equation during withdrawal depends on the maturation flux (k_maturation·N_silent) rather than the absolute count of mature synapses. This has important implications:

The computational framework generates several testable predictions:

The model identifies distinct intervention windows (Figure 4B):

The model employs deterministic dynamics, simplified discrete states, and heuristic plasticity weights. The parameter sensitivity analysis demonstrates robustness across ±50% variations, and alternative initial conditions (±10%–20% juvenile at baseline) did not alter qualitative conclusions. Future refinements should incorporate stochastic dynamics, spatial organization, and systematic parameter estimation from experimental data.