Computational models of dopamine release measured by fast scan cyclic voltammetry in vivo

Abstract Dopamine neurotransmission in the striatum is central to many normal and disease functions. Ventral midbrain dopamine neurons exhibit ongoing tonic firing that produces low extrasynaptic levels of dopamine below the detection of conventional extrasynaptic cyclic voltammetry (∼10–20 nanomolar), with superimposed bursts that can saturate the dopamine uptake transporter and produce transient micromolar concentrations. The bursts are known to lead to marked presynaptic plasticity via multiple mechanisms, but analysis methods for these kinetic parameters are limited. To provide a deeper understanding of the mechanics of the modulation of dopamine neurotransmission by physiological, genetic, and pharmacological means, we present three computational models of dopamine release with different levels of spatiotemporal complexity to analyze in vivo fast-scan cyclic voltammetry recordings from the dorsal striatum of mice. The models accurately fit to cyclic voltammetry data and provide estimates of presynaptic dopamine facilitation/depression kinetics and dopamine transporter reuptake kinetics, and we used the models to analyze the role of synuclein proteins in neurotransmission. The models’ results support recent findings linking the presynaptic protein α-synuclein to the short-term facilitation and long-term depression of dopamine release, as well as reveal a new role for β-synuclein and/or γ-synuclein in the long-term regulation of dopamine reuptake.


Introduction
Dopamine (DA) is a neurotransmitter that plays important roles in learning and memory, drug and alcohol addiction, attention deficit disorder, schizophrenia, and neurodegenerative disorders including Parkinson's disease (1). Electrochemical recording methods, particularly fast-scan cyclic voltammetry (FSCV) and constant potential amperometry, are used to record and analyze evoked DA release from the axons of dopamine neurons in rodent brains in vivo (2,3), ex vivo in brain slice preparations (4,5), and in cultured neurons (6,7). The development of computational models to quantitatively describe FSCV and amperometry recordings of evoked DA release (8)(9)(10)(11)(12)(13) has enabled the detailed analysis of DA release and reuptake kinetics under different stimulation paradigms both ex vivo and in vivo. However, the plasticity of DA neurotransmission is far more complex in vivo due to multiple factors that are less important in culture or slice preparations, including modulatory inputs from neurotransmitters such as acetylcholine and γ-aminobutyric acid (GABA) and multiple presynaptic regulatory proteins such as the synapsins and synucleins. Importantly, ventral midbrain DA neurons are tonic pacemakers, but their cell bodies are absent in slice preparations.
In vivo, tonic firing produces minimal DA levels in the striatum and is interspersed with burst-firing activity that can cause DA release to saturate the dopamine transporter (DAT) reuptake system, thus providing a nonlinear increase in extracellular DA during bursts (14). As reported in a study of spontaneous DA release (9), the magnitude and duration of burst-firing-evoked DA release in vivo are subject to multiple depressing and facilitatory influences that remain mostly uncharacterized. As such, deciphering the kinetics of DA release during burst firing is central for elucidating the interactions of DA in vivo and providing insights into the mechanisms responsible for DA release and regulation in normal development and function, exposure to drugs (ranging from antipsychotics to recreational/abused drugs), and many disease states. We have created three computational models of burst-firing striatal DA release in vivo, each of which closely fits the data and can be used to derive biologically interpretable parameters such as estimates of facilitation/depression kinetics, amount of DA release, and DAT reuptake kinetics. These models account for multiple biological, electrochemical, and experimental factors that affect the FSCV trace and can be adapted for other approaches to detect extracellular DA (15) and for analyzing release kinetics of other neurotransmitters.
The models simulate the interactions between extracellular DA released into striatal tissue, the carbon-fiber electrode that detects DA, and an intermediate area of damaged tissue known as the "dead space," formed when inserting the electrode into the tissue (16) (Fig. 1A). When an electrical stimulus is applied, DA released by vesicular exocytosis diffuses through the extracellular space and dead space until it reaches the carbon-fiber electrode. In contrast to amperometry electrodes that "consume" DA molecules through oxidation, FSCV electrodes act as reflective surfaces that do not consume DA molecules and cause them to diffuse back into the tissue through the dead space (5). As a result, the concentration of extracellular DA in the striatum with FSCV experiments is mostly regulated by DAT uptake. The concentration of DA measured at the electrode is further altered by electrochemical adsorption (17), which causes the signal baseline to increase when multiple bursts of stimuli are triggered in quick succession. Our models can also be used to simulate DA release and diffusion without the recording "artifacts" introduced by the electrode, dead space, and adsorption (Fig. S1).
While the three models described here differ in the implementation of release and diffusion mechanics, their functional predictions are very similar. The Simple Uniform Release (SUR) model ( Fig. 1B) is one-dimensional and computes DA release as a uniform distribution in the striatum that diffuses toward the carbon-fiber electrode without taking into account physical characteristics of the striatum. The Spatiotemporal Uniform Release (STUR) model ( Fig. 1C) is two-dimensional and incorporates spatiotemporal diffusion of DA in the striatum using cylindrical coordinates and a uniform release distribution similar to the SUR model. The Spatiotemporal Discrete Release (STDR) model (Fig. 1D) is similar to the STUR model but uses discrete DA release sites positioned throughout the striatum. The mathematical equations for these models are described in the Materials and Methods section.
Our models contribute to the current state-of-the-art research in four key aspects. First, our models describe plasticity kinetics of evoked DA release for in vivo FSCV traces from anesthetized mice, in contrast to a similar study primarily focused on awake mice (9). Second, the spatiotemporal models provide mathematical equations for discrete release sites and the dead space, whereas prior studies (10, 11) mention these features without equations. Third, the creation of three computational models with different release and diffusion mechanics gives researchers the flexibility to choose the model that best captures their specific experimental conditions. Finally, we provide a statistical inference algorithm known as automatic differentiation variational inference (ADVI) (19,20) for researchers to fit the adjustable parameters (with associated uncertainty estimates) in the models based on biological constraints.

Results
We fit the computational models to in vivo recordings of dorsal striatum DA release evoked by multiple burst stimuli in the ventral midbrain (21). The complete list of model parameter estimates is provided in the Supplementary Information (Tables S1,  S2, S3). Briefly, the experimental procedure was comprised of 6 "Single Burst" and "Repeated Burst" stimulation protocols (each pass being termed a "sweep"), with a 2-min recovery time between protocols and a 6-min recovery time between sweeps ( Fig. 2A). The Single Burst protocol was a train of 30 pulses at 50 Hz, and the Repeated Burst protocol was a series of six Single Bursts with a 5 s interstimulus period between two bursts. The relationship between number of pulses and pulse frequency to DA release is linear with these stimulus parameters (21), as well as the relationship between stimulation current and DA release (Fig. 2B).

Three kinetic components are sufficient to describe DA release kinetics in wildtype mice in vivo
To analyze the kinetics of DA release under normal conditions, the models were fit to FSCV data from wildtype (WT) mice using the Single Burst and Repeated Burst protocols from Sweep 1 (Fig. 3) and Sweep 6 ( Fig. 4) of the experimental procedure. Consistent with conclusions from a study of spontaneous DA release kinetics in behaving mice (9), the complex dynamics of the data require three kinetic components to capture changes in DA release over time: short-term facilitation, short-term depression, and longterm depression. Each kinetic is described in the models by a "plasticity factor" p j that controls the facilitative/depressive effect's magnitude and a time constant τ j that controls the kinetic's duration (see Materials and Methods for more information). The kinetic parameter estimates (Table 1) for the Single Burst and Repeated Burst protocols indicate that in WT mice, the plasticity factor for short-term facilitation is 3.5-fold larger than short-term depression, while the time constant of short-term depression is ∼1.7-fold to 2-fold longer than short-term facilitation.
The long-term depression kinetic contributes to an ∼50% decrease in evoked DA release observed in WT mice between Sweeps 1 and 6. While the plasticity factor of long-term depression is weaker than both short-term kinetics (Table 1), the time constant of 15 min is ∼60-fold longer than the closest time constant for short-term kinetics, consistent with prior literature (9,22). The release parameter estimates (Table 1) indicate that DAT activity V m and DA release DA P decrease between Sweeps 1 and 6 by 33.3% and ∼25-30%, respectively. Further experiments might investigate the biological mechanisms that relate to long-term plasticity of release and reuptake that could account for these effects, such as synaptic exhaustion at vesicular release sites and DAT endocytosis.

Synuclein proteins regulate DA release kinetic components in vivo
In order to determine how the kinetic components can be altered in vivo, we analyzed the impact of the synuclein family of neuronal proteins on DA release. There are three genetically distinct synuclein proteins: α-synuclein (α-Syn), β-synuclein (β-Syn), and γ-synuclein (γ-Syn) (23). α-Syn and β-Syn are expressed in brain regions including the striatum, thalamus, hippocampus, neocortex, and cerebellum (24), and γ-Syn is expressed in the peripheral nervous system as well as the brain (25). α-Syn has been associated with disease-related protein aggregates in Parkinson's disease and other neurodegenerative disorders, β-Syn has been shown to inhibit α-Syn aggregation, and abnormal levels of γ-Syn have been linked to multiple cancers (26). Although the normal physiological functions of the synucleins are still under study, they regulate evoked DA release both in the striatal slice (27) and in vivo (21), likely by modulating steps in synaptic vesicle recycling and exocytosis (28,29).
To analyze the interactions of synucleins on DA release in vivo and its effects on the kinetics of DA release, the models were fit to FSCV data from two lines of mutant mice: α-Syn knockout (KO) mice (Figs. 5 and 6) that are deficient in α-Syn, and synuclein triple knockout (Syn TKO) mice (Figs. 7 and 8) that are deficient in α-Syn, β-Syn, and γ-Syn. Like WT mice, the data were taken from the Single Burst and Repeated Burst protocols for Sweep 1 and Sweep 6. Comparing the models' kinetic parameter estimates for α-Syn KO mice (Table 1) revealed that the plasticity factor of short-term facilitation decreased by ∼40-50%, indicating that α-Syn expression enhances short-term facilitation. Additionally, the time constant of short-term depression increased by ∼2.5-fold, indicating that α-Syn expression slows down the kinetics of short-term depression. Finally, the plasticity factor of long-term depression is absent, implying that long-term depression is fully dependent on α-Syn expression. Similar changes in facilitation/depression constants were observed in Syn TKO mice, indicating that the effect was not related to β-Syn or γ-Syn expression. These results reinforce conclusions by Somayaji et al. (21) that α-Syn controls the short-term facilitation and long-term depression of DA release.
Similar to WT mice, the models' release parameter estimates (Table 1) for α-Syn KO mice indicate a 33.3% decrease in DAT activity and ∼30-35% decrease in DA release between Sweeps 1 and 6, which accounts for the smaller ∼25% decrease in evoked DA release observed in the data between the sweeps in the absence of any long-term depression kinetic. In contrast, the models' estimates for the Syn TKO mice indicate a relatively small increase in both DAT activity and DA release between Sweeps 1 and 6 by 16.7% and ∼5-10%, respectively (Table 1). This suggests a previously unknown role of β-Syn and/or γ-Syn in the long-term regulation of DAT uptake and vesicular DA release and may in part underlie the observation that the Syn TKO line has a greater effect than the α-Syn KO line on exocytosis of peptides from chromaffin cells and cultured hippocampal neurons (30).

Discussion
In this paper, we presented three computational models of DA release to fit and analyze FSCV recordings of evoked DA release in the dorsal striatum of anesthetized mice in vivo. Using these models, we quantified biochemical mechanisms that control DA neurotransmission, including the kinetics of DA release, DAT uptake, and DA diffusion. Similar to a previous analysis of DA

Electrode Dead Space
[DA] η [DA] η P Fig. 1. Schematic representations of computational models. A, Coronal slice of the mouse brain (adopted from the Allen Reference Atlas -Mouse Brain (18)). During in vivo FSCV experiments, a bipolar electrode in the ventral midbrain (not shown) is used to electrically stimulate the cell bodies of dopaminergic neurons, causing DA release via vesicular exocytosis from dopaminergic axons in the dorsal striatum (blue). The DA molecules diffuse to the carbon-fiber electrode (black) through striatal tissue and an area of damaged tissue known as the "dead space" (orange) that is formed when the electrode is inserted into the brain. B, Simple Uniform Release Model. DA release sites are modeled as a continuous release source within the striatum using temporal diffusion in one dimension without calculating the spatial characteristics of the striatum. C, Spatiotemporal Uniform Release Model, which is similar to the Simple Uniform Release model, except that the spatial characteristics of the dorsal striatum are modeled as a cylinder using isotropic diffusion in two dimensions. D, Spatiotemporal Discrete Release Model, which is similar to the Spatiotemporal Uniform Release model, with the addition of release sites modeled as discrete point sources positioned throughout the striatum.
neurotransmission in awake mice (9), we found that three kinetic components effectively describe modulation of DA release over time: short-term facilitation, short-term depression, and long-term depression. We also analyzed the interactions of the synuclein family of proteins with DA release and mathematically quantified their effect. Our results support previous findings that α-Syn controls the short-term facilitation and long-term depression of DA release. Furthermore, we identified a new role for β-Syn and/or γ-Syn in the long-term regulation of DAT uptake and DA release. The computational models we present have only been tested for in vivo FSCV recordings but could be adapted to ex vivo recordings in the striatal slice. The rapid recovery time of evoked DA release in vivo of < 10 ms (21) allows the electrical stimulation of midbrain neurons at 50 Hz, providing the means to study short-term synaptic plasticity. In contrast, the recovery time between electrical stimuli in slice is ∼1 min (31), which is due in part to the lack of tonic DA neuron firing and the lack of other neuronal activity in slice except for the cholinergic interneuron (32). As such, single electrical pulses cannot be recorded in the striatal brain slice without an extended recovery time and rapid electrical pulse trains produce a negligible increase in DA release.
Given that the three computational models produce similar fits to the data, they can be used without fundamentally altering the model's estimates for the DA release kinetics or DAT uptake. For fitting data with multiple burst stimuli such as the Repeated Burst protocol, users may achieve slightly improved computational performance using the Simple Uniform Release Model. However, it may be beneficial to use the spatiotemporal models if the FSCV experiments are performed in conjunction with other electrochemical or optical procedures to elucidate physical characteristics in the striatum. For example, it is possible to use amperometry to estimate the size of the dead space around the electrode, and the current dead space radius of 3 μm specified for Figs. 3-8 was derived from an amperometric study (10). As the dead space may vary from experiment to experiment which will alter the DA kinetics measured by the carbon-fiber electrode (Fig. 9), using amperometry and FSCV together can enable the dead space dimensions to be computed from amperometry data and used in the spatiotemporal models to fit FSCV data.
Additionally, nonelectrochemical DA recording methods such as optical reporters or false fluorescent neurotransmitters (FFNs) can be used to directly visualize the distribution of active and silent release sites in the striatum, which would make the Spatiotemporal Discrete Release Model particularly beneficial compared to other models. Currently, the distance between two adjacent release sites in this model for the fits presented in Figs. 3-8 is set to 6 μm based on numbers derived from FFN imaging reported by Pereira et al. (33) but assumes that the release sites are spaced at equal intervals. In theory, if the FSCV experiments are performed in conjunction with FFN or optical reporter experiments, the position and distance of release sites can be adjusted to any location according to the experiment. While our analysis of FSCV provides genuine measures of DA concentration, in the future, it may be possible to adapt the models and describe  kinetics of the change in fluorescence (ΔF/F ) data from fluorescent reporters, as many labs have access to fiber photometry. The cylindrical diffusion of DA implemented in the spatiotemporal models may be particularly useful for modeling DA release under the influence of DAT-inhibiting drugs including amphetamine, cocaine (34), and nomifensine (35). The diameter of the cylinder is set to 100 μm for the fits presented in Figs. 3-8, less than the actual mouse striatum diameter of 2-3 mm (36). Under typical circumstances, distant release sites do not contribute to the final FSCV trace because DAT will remove distally released DA before it diffuses to the electrode. However, if DAT is at low levels or inhibited through pharmacological methods, the computed diameter of the striatum can be adjusted to capture DA diffusion in the extracellular space from distant release sites and account for the increase in DA concentration at the electrode.
With multiple adjustable parameters in the computational models that can be used to alter the outputs, it is theoretically possible to generate multiple sets of "best-fit" solutions that provide close fits to the FSCV data. In practice, however, the range of solutions is limited by the biological constraints applied to each parameter. For example, the affinity/binding constant of DAT (K m ) can theoretically be any nonnegative value, but a consensus in the research literature is for K m = 0.2 μM under normal conditions in mice (8), so fits that deviate from this standard without a valid biological reason can be ignored. To assist with providing the best fits with our computational models, we implemented a statistical inference algorithm known as ADVI (19,20) for curve fitting that can limit the ranges for each parameter to effectively factor in the biological constraints.

Computational models
The three models presented in this paper were implemented in Python, and the detailed discretization schemes used in the implementation can be found in the Supplementary Information. We used a combination of ordinary differential equations (ODEs), partial differential equations (PDEs), and time-series functions to implement these models.

Stimulation pattern
The function S is used to compute the electrical bursts applied from the stimulating bipolar electrode in FSCV experiments with the Heaviside theta function θ: In this step function, S(t) = 1 denotes a time where an electrical burst is present, and S(t) = 0 denotes a time where an electrical burst is absent. The variable t i denotes the start time for each burst, NP is the number of pulses in each burst, and f is the stimulus frequency (in Hz). Due to the increased lag time introduced by spatial diffusion in the dead space, t i is modeled earlier in the STUR and STDR models compared to the SUR model.

Plasticity kinetics
The kinetics of stimulation-dependent DA release A is computed as the product of individual kinetic components H j , revised from Montague et al. (9): During stimulus periods when S(t) = 1, a multiplicative plasticity factor p j is applied to each kinetic component; positive values for p j produce facilitation and negative values for p j produce depression. During interstimulus periods when S(t) = 0, the kinetic components decay toward the equilibrium value of 1 based on the time constant τ j until the next electrical burst is triggered. An example of the individual kinetic components for short-term facilitation, short-term depression, and long-term depression is shown in Fig. S2.

DA released into the striatum
Using the plasticity kinetics A and stimulation pattern S defined above, the DA released in the striatum [DA] S is computed as the net difference between the concentration of DA released by electrical stimulation into the striatum and removed through DAT uptake. Each model implements [DA] S differently based on the different temporal and spatial characteristics.
In the Simple Uniform Release Model, [DA] S is implemented as an ODE that computes the change in concentration over time: DA P is the amount of DA released per milliamp (mA), and I is the electrical stimulus current in mA, such that DA P I is the concentration of DA released per pulse. The linear relationship established between DA release and stimulation current with DA P I can be found in FSCV experiments that apply electrical stimulation between 0.3 and 0.6 mA (Fig. 2B). V m and K m are the maximal velocity and affinity constant of DAT uptake (37), modeled using firstorder Michaelis-Menten kinetics (38).
In the Spatiotemporal Uniform Release Model, [DA] S is implemented as a PDE that simulates spatiotemporal diffusion and assigns a physical dimension to the dead space in brain tissue: The striatum is modeled as a cylinder with uniformly distributed release sites, and the carbon-fiber electrode and dead space are positioned in the center of the striatum. The diffusion of DA in the striatum is simulated as isotropic diffusion using radial coordinates as described by John Crank (39), where R signifies the distance along the radius of the cylinder and D = 240 μm 2 /s is the diffusion coefficient of DA in the brain accounting for tortuosity (40). The dead space η is computed using the Heaviside theta function θ, where R L is the radius of the cylinder and R D is the radius of the dead space. The area covered by the dead space does not include DA release and DAT uptake due to the damaged tissue (10). In the Spatiotemporal Discrete Release Model, [DA] S is implemented nearly identically to the Spatiotemporal Uniform Release Model, except for the inclusion of an additional variable P to simulate spatially discrete release sites instead of a continuous release area: P uses the Dirac delta function δ to simulate multiple "point" sources of release in the cylinder, with R i specifying the location of each point along the radius, and each point source is a concentric cylinder due to the radial diffusion. Depending on the distribution of release points, the estimates for DA P can be higher compared to the Simple Uniform Release and Spatiotemporal Uniform Release models (see Table 1); sparser distributions lead to higher release values from each individual point source.

DA measured at the electrode
Finally, the measurement of DA at the electrode is implemented as a coupled ODE system to model the temporal response and electrochemical adsorption that occur with the carbon-fiber electrode in FSCV experiments: [DA] E computes the DA diffusing from the extracellular space based on the value computed for [DA] S minus the DA leaving the electrode due to the reflective surface in FSCV experiments, plus the residual DA that occurs due to adsorption on the electrode. k S is the rate transfer of DA moving from the striatum toward the electrode, k E is the rate transfer of DA moving away from the electrode, and k Γ = 1 is the rate transfer of DA adsorbing  to the electrode. The adsorption Γ DA is modeled as the concentration of DA that adsorbs to the electrode minus the concentration of DA that desorbs from the electrode, and it is a modified version of an equation presented by Bath et al. (17). The adsorption kinetic k ads 1 controls the amount of DA that adheres to the electrode, and the desorption kinetics k ads 2 and k ads 3 control the amount of residual DA on the electrode that falls off. The full derivation for Equation (11) is provided in the Supplementary Information.

Curve fitting
We implemented a statistical inference algorithm for curve fitting known as ADVI (19,20) to assist with optimizing the parameters in our computational models and computing the best fit to the data. ADVI differs from conventional curve fitting algorithms used for parameter estimation in a few key aspects. First, ADVI can define constraints for each parameter estimate within a certain range of values, which is especially important for the models as it ensures biologically reasonable values that can still closely fit the data. Second, ADVI returns mean and standard deviation estimates for each parameter upon optimization rather than a singular point value; the standard deviation provides a measure of "uncertainty" for each parameter estimate and shows how much it could potentially vary while still explaining the data well. Third, ADVI is more efficient as it is a gradient-based optimization approach that can compute the optimization step in one pass for each iteration, in contrast to traditional nonconvex optimization approaches which require repeated numerical integration or Monte Carlo sampling.
The algorithm works by transforming the parameters from the original latent variable space to the unconstrained latent variable space. It then uses both the unconstrained and original parameters to compute an objective called the evidence lower bound (ELBO), comprised of the expected log joint density, the log Jacobian of the transformation, and the entropy of the unconstrained parameters (which are drawn from a Gaussian distribution) (19). The gradients of the unconstrained parameters with respect to the ELBO are used to optimize the original parameters, and the algorithm iterates until the parameters converge. Because the gradients can be computed efficiently using automatic differentiation, variational inference is typically faster than other statistical inference methods like Markov Chain Monte Carlo (MCMC) (41).
Applying ADVI requires adapting the ODE model into a probabilistic model to compute the expected log joint density for the ELBO, comprised of the log-prior and log-likelihood. To calculate the log-prior, the free parameters in the ODE model (V m , [DA] P , p j , τ j , k S , k E , k ads 1 , k ads 2 , and k ads 3 ) are drawn from independent Gaussian distributions. For the log-likelihood, the ODE's solution at a given state Ĉ t is computed by integrating [DA] E over time, shown in Equation (12). Ĉ t is then assigned as the mean of a Gaussian distribution to compute the log-likelihood (with the standard deviation v set to 0.05), and the corresponding state in the real data C t is drawn from that distribution, as in Equation (13):

In vivo recordings
To validate the computational models, we used previously published in vivo FSCV traces of DA release from WT and transgenic mice (21). Briefly, male or female mice were anesthetized with isoflurane (induction 2.5%, maintenance 0.8-1.4% in O 2 , 0.35 L/min), placed on a heating pad, and head-fixed on a stereotaxic frame. A craniotomy was performed to insert a 22G bipolar stimulating electrode (P1 Technologies) in the ventral midbrain to trigger electrical pulse trains and a carbon-fiber microelectrode (5 μm diameter, ∼150 μm length) in the dorsal striatum to measure evoked DA release; detailed stereotaxic coordinates and procedures can be found in the publication. The depth of the stimulating electrode was adjusted between 4 and 4.5 mm for maximal DA release. To conduct FSCV, a triangular voltage waveform (−450 to +800 mV at 294 mV/ms versus Ag/AgCl reference electrode) was applied to the carbon-fiber electrode at a 10 Hz sampling rate (or every 100 ms). The signals measured at the electrode were converted from analog to digital signals using a 16-bit data acquisition interface and subsequently recorded using IGOR Pro software. The electrical pulse trains were delivered to the stimulating electrode at a constant current of 400 μA using a stimulus isolator and a pulse generator. The stimulation protocols used in the experiments are explained in the Results. The carbon-fiber microelectrodes were calibrated in artificial cerebrospinal fluid using known concentrations of DA.

Supplementary Material
Supplementary material is available at PNAS Nexus online.

Funding
This research was supported in part by NSF GRFP grant DGE 2036197 (N.S.), European Research Council Synergy grant ASTRA  Table S2).