A model for NMDA-receptor-mediated synaptic currents generating persistent activity proposed by Wang and Brunel [1–3] has been widely adopted in computational neuroscience, both for spiking-neuron and mean-field models [2, 4]. The model describes NMDA synaptic dynamics by a phenomenological two-dimensional nonlinear ODE system for the gating variable $S(t)$ evolving according to $$ \frac{dS_j}{dt} = -\frac{S_j}{\tau_d} + \alpha x_j (1 - S_j) \\ \frac{dx_j}{dt} = - \frac{x_j}{\tau_r} + \sum_k \delta(t - t_j^k) $$ where indices $j$ mark presynaptic neurons and $k$ spike times. Due to the nonlinearity, the pre-synaptic gating variables of a post-synaptic neuron cannot be simulated in aggregated form. Numerically efficient solutions are only practical for networks with simple connectivity structure and identical, short delays as in a Brian2 implementation [5] of the network studied in [2]. Noting that $S_j(t)$ only depends on the spike history of the presynaptic neuron, we show that the NMDA dynamics can be approximated by an exponential decay between spikes and a history-dependent jump upon spikes. Let $t_{\text{ls}}$ be the time of the previous spike of neuron $j$ and $t$ of the current spike. We then set $$ S_j(t^+) = k_0 + k'_1 S_j(t^-) $$ where $$ \begin{align} S_j(t^-) &= S_j(t_{\text{ls}}) e^{- \frac{t^- - t_{\text{ls}}}{\tau_d}} \\ k_0 &= (\alpha \tau_r)^{\frac{\tau_r}{\tau_d}} \gamma \left[1 - \frac{\tau_r}{\tau_d}, \alpha \tau_r \right] \\ k'_1 &= e^{-\alpha \tau_r} \end{align} $$ with the lower incomplete gamma function $\gamma(a, z)$. This model is scheduled for inclusion in NEST 3.8 as `iaf_bw_2001`. As shown in the figure, this model approximates the exact solution well, including the decision-making behavior reported by Brunel and Wang. The computational benefit of the approximate model grows rapidly with the size of the network, and for a fully connected network of 16,000 neurons run with 12 threads, the approximation is about 1,060 times faster than the exact model in NEST. With the same setup, the approximation is also roughly 3 times faster than the previously mentioned efficient solution for fully connected networks with identical, short delays, implemented in Brian2. Exploiting the flexibility and performance gained through the approximation, we investigate the dynamics of a binary decision-making network with sparse connectivity and randomized delays.