This repository (repo) contains the codes for a manuscript titled Sequential Transitions of Male Sexual Behaviors Driven by Dual Acetylcholine-Dopamine Dynamics submitted to Cell press. For reproducibility, the codes and data in this repo can reproduce all the reported results in our paper.
(NT) in the brain is crucial for understanding of the precise neural basis. Here, we reveal what kind of neurons and receptors are important for the specific NT release by computational modeling based on the real NT and neuronal activity patterns in the Nucleaus accumbuns. As we observed similar dynamics in NT release and neuronal activity through GRAB and GCaMP fiber photometry imaging (GRABDA2m and GCaMP6s of dopaminergic axons in the vsNAc, GRABACh3.0 and GCaMP6s of cholinergic neurons in the vsNAc), we hypothesized that (1) the neuronal activity of a ensemble of neurons reflects the amount of NT released; (2) the actual GCaMP recordings can predict the release probability of the NT. Based on these hypotheses, we developed a computational model to estimate the cellular environment in which neurons and receptors are involved in generating dual ACh-DA rhythms during intromission."
Refering to the code/run
, execute python file such as experiment.py
in line 22.
-
code/resources
: definition of the hps space-
exp01.yaml
: To define the hps space -
trial.csv
: To define$Activation$ pattern ($ChAT_{sim}^{vsNAc}$ /$DA_{sim}^{VTA→vsNAc}$ = on/on) and the receptor distribution (all 4 receptors are distributed) to visualize hps -
trial0.csv
: To define$Activation$ pattern and the receptor distribution for statistics
-
-
code/src
: basic codes-
cells.py
: To define each element involved in the neurotransmitter (NT) release. -
data.py
: For loading data. -
simulator.py
: For simulating$NT_{sim}$ release. -
utils.py
: To make the code easier to follow, all basic codes that are not directly relevant to the analysis logics are inside utils.py.
-
-
code
: source codes for running analysis-
config.yaml
: To define the data path and basic information of$NT_{real}$ and$NT_{sim}$ -
experiment.py
: For running the simulation experiment. -
hps_visualization.py
: For visualization of$Diff(NT_{sim}, NT_{real})$ and hps. -
run
: To perform modeling.- It takes 3 hours for the simulation for one setting. You can change the search time in
exp01.yaml
.
hps_config: max_time: 10800 #(sec)
- You can define the required auguments (in
code/resources
) inexperiment.py
by argparse inrun
.
# for hps search for statistics using trial0.csv python experiment.py --setting resources/trial0.csv
- It takes 3 hours for the simulation for one setting. You can change the search time in
-
run_trial.py
: For trial run to check if the simulator works. -
run_hps.py
: For executing single hyperparameter search (hps). -
run_eval.py
: For evaluating specific simulation setting.
-
-
data
: typical rhythmic signals of in vivo GRABACh3.0 and rGRABDA2m during intromission ($NT_{real}$ ) -
results
: analysis results including hps and the simulated signal ($NT_{sim}$ )-
eval
: Comparison of$NT_{sim}$ and$NT_{real}$ -
hps
:$Diff(NT_{sim}, NT_{real})$ and hps.
-
- Required python libraries is listed under
code/requirements.txt
.
- Distributed under the MIT License. See
LICENSE
for more information.
-
In this model, we set neuronal assembly (
$ChAT_{sim}^{vsNAc}$ asclass AChNeuron
and$DA_{sim}^{VTA→vsNAc}$ asclass DopamineAxon
in thecode/src/cells.py
), the internal activity of neurons (release probability measurement,$RPM$ asrelease_prob_idx
in thecode/src/cells.py
), NT release ($ACh_{sim}$ and$DA_{sim}$ ) based on the$RPM$ (calculated bydef release
in eachclass AChNeuron
andclass DopamineAxon
in thecode/src/cells.py
), receptors expressed on neurons (D1R, D2R, and nAChR), and the extracellular environment in which the NT can be released. -
To mimic the change of the amount of the extracellular NT over time, in this model, we assumed the three steps of cycle (as
def step
in eachclass Environment
in thecode/src/simulator.py
) :- (i) the activation of target neurons by the external stimuli (
$Activation$ asdef Activation
in eachclass AChNeuron
andclass DopamineAxon
in thecode/src/cells.py
andclass Activator
in thecode/src/utils.py
); - (ii) NT released from the target neurons (as
def release
in eachclass AChNeuron
andclass DopamineAxon
in thecode/src/cells.py
); - (iii) the effect of the receptor by binding NT to the receptor on the target neurons (as
def receptor_nAChR
,def receptor_d1
, anddef receptor_d2
inclass AChNeuron
andclass DopamineAxon
in thecode/src/cells.py
).
- (i) the activation of target neurons by the external stimuli (
-
These cycles were repeated a certain number of times and, during this period, we recorded the
$RPM$ and amount of$NT_{sim}$ released from the neurons which was determined based on the$RPM$ value. Through this model, we attempted to estimate the appropriate cellular environment, including the type of receptors or neurons involved, by minimising$Diff(NT_{sim}, NT_{real})$ , the difference between the actual ACh-DA dynamics in vivo (NT recorded by fiber photometry:$NT_{real}$ ) and model-simulated ACh-DA dynamics (NT generated by the model:$NT_{sim}$ ) (Figures 4A, 4B, and S4A). -
Hypothetical neurons and receptors used in this simulation were identified based on ex vivo imaging (Figure 3), the in situ hybridisation results (Figure S3), and previous reports
$^1$ . As neurons and receptors, we placed D1R and D2R in$ChAT_{sim}^{vsNAc}$ neurons (asdef receptor_d1
anddef receptor_d2
inclass AChNeuron
in thecode/src/cells.py
, respectively), as well as D2R and nAChR in$DA_{sim}^{VTA→vsNAc}$ axon terminals (asdef receptor_d2
anddef receptor_nAChR
inclass DopamineAxon
in thecode/src/cells.py
, respectively).
- To mimic the cellular environment during neurons are activated by the external stimuli, in this modeling, we assumed that the intensity of neuronal activity was changed by the external stimuli or the effect of the receptor expressed on the neurons.
$Activation$ includes the spontaneous activity of the target neurons and stimulation of neurons by substances other than ACh and DA. Here, we defined that the$Activation$ continuously adds a certain value to the target neuron (class Activator
in thecode/src/utils.py
). Thus, the value added as an external stimulus to the$RPM$ of the target neurons every cycle was constant.
- To mimic the manner of NT release, we defined that the volume of
$NT_{sim}$ released to the extracellular environment from the neurons was determined by the$RPM$ value of the neurons. In this model, the volume of$NT_{sim}$ released is equivalent to the$RPM$ value if$RPM$ is positive and 0 if$RPM$ is negative (Figure 4B) (def relu
in thecode/src/utils.py
). Additionally, we defined the diffusion function of DA to mimic the in vivo extracellular DA diffusion mechanism (def da_split
anddef da_diffusion
inclass Environment
in thecode/src/simulator.py
)$^2$ . The time course and intensity of the DA diffusion were investigated using this model (Figures S4B-S4D).
-
Next, we considered a strategy to mimic the behavior of nAChR, D1R, and D2R. As reported in previous studies, DA receptors are GPCRs that are thought to respond more slowly than nAChR, the ionotropic receptor
$^{3–7}$ . Incorporating these insights, we defined the characteristics of each receptor using two variables:$RE$ (receptor efficacy), the effect intensity which changes$RPM$ ;$RS$ (receptor speed), the latency of the effects occurs in$RPM$ after$NT_{sim}$ binds to the receptor. -
In this experiment, to simplify the setting, we set the nAChR related variables as constant values (
$nAChR's RE =1$ asnAChR_efficacy
, and the$RS =1$ ), and estimated the characteristics of D1R and D2R by exploring their$RE$ and$RS$ .$RE$ of D1R and D2R was searched over a range of positive values, and their$RS$ was searched over a range of natural numbers.- For example, when
$ACh_{sim}$ binds to the nAChR expressed on the$DA_{sim}^{VTA→vsNAc}$ axon terminals, there is an increase in the$RPM$ of the$DA_{sim}^{VTA→vsNAc}$ axon terminal$ACh_{sim} × 1$ ($nAChR's RE =1$ ) one step later ($nAChR's RS =1$ ). Since the D1R receptor is an excitatory GPCR, we defined its operation as follows: when D1R receives$DA_{sim}$ ,$D1R's RS$ (asdelay_d1
in thecode/src/cells.py
) step later, the$RPM$ is increased to$DA_{sim} × D1R's RE$ (asd1_efficacy
in thecode/src/cells.py
). - By contrast, we defined operation of the D2R receptor, an inhibitory GPCR, receiving as
$D2R's RS$ (asdelay_d2
in thecode/src/cells.py
) step later,$RPM$ is decreased by a factor of$DA_{sim} × D2R's RE$ (asd2_efficacy
in thecode/src/cells.py
). With these receptors' settings, we conducted computational modeling to estimate the$RE$ and$RS$ values.
- For example, when
-
To estimate receptor efficacy (
$RE$ ) and receptor speed ($RS$ ) of DA receptors, we explored$RE$ and$RS$ values which minimise$Diff(NT_{sim}, NT_{real})$ , the difference between$NT_{sim}$ and NTreal using Bayesian optimisation (Figure 4B) (code/run_hps.py
). -
We used Optuna
$^8$ , an automatic hyperparameter optimisation software framework, to determine the best values for$RE$ and$RS$ using$Diff(NT_{sim}, NT_{real})$ as an indicator. We used the following values as a fixed value throughout exploration. The value added by$Activation$ is 1. The number of trials conducted to estimate$RE$ and$RS$ of D1R and D2R was approximately 1,000 to 10,000, depending on the experimental design. In each experimental design, out of these trials, the one representing the smallest$Diff(NT_{sim}, NT_{real})$ , the condition where$NT_{sim}$ and NTreal are considered as suitable values, was used to define$RE$ and$RS$ of D1R and D2R. To express the effect speed of each receptor in seconds, the RS of each receptor and DA diffusion step were divided by 200, the sampling rate (Figures 4E and S4D).
- Ai Miyasaka
- University of Tsukuba
- Naoki Nonaka
- RIKEN
- Takeshi Kanda
- Nara Medical University
- Yuka Terakoshi
- University of Tsukuba
- Yoan Cherasse
- University of Tsukuba
- Yukiko Ishikawa
- University of Tsukuba
- Yulong Li
- Peking University
- Hotaka Takizawa
- University of Tsukuba
- Jun Seita
- RIKEN
- Masashi Yanagisawa
- University of Tsukuba
- Katsuyasu Sakurai
- University of Tsukuba
- Takeshi Sakurai
- University of Tsukuba
- Qinghua Liu
- National Institute of Biological Sciences, Beijing
corresponding_contributors:
- Qinghua Liu ([email protected])
- Takeshi Sakurai ([email protected])
- Katsuyasu Sakurai ([email protected])