Year: 2021

Abstract:
Brain-wide measurements of neural activity combined with detailed measurements of anatomical connectivity of the C. elegans nervous system in principle allow for the development of detailed mechanistic computational models. However, there are several challenges. We often do not have direct experimental access to important modeling details such as single-neuron dynamics and the signs and strengths of the synaptic connectivity. Further, neural activity can only be measured in a subset of neurons, often indirectly via calcium imaging, and significant trial-totrial variability has been observed. To overcome these challenges, we introduce a connectome-constrained latent variable model (CC-LVM) of the unobserved voltage dynamics of the entire C. elegans nervous system and the observed calcium signals. We used the framework of variational autoencoders to fit parameters of the mechanistic simulation constituting the generative model of the LVM to calcium imaging observations. A variational approximate posterior distribution over latent voltage traces for all neurons is efficiently inferred using an inference network, and constrained by a prior distribution given by the biophysical simulation of neural dynamics. When applied to a recent dataset, we find that connectomic constraints enable our LVM to predict the activity of neurons whose activity were withheld. We explored models with different degrees of biophysical detail, and found that models with the most realistic conductance-based synapses provide markedly better predictions than current-based synapses for this system.

text: Brain-wide measurements of neural activity combined with detailed measurements of anatomical connectivity of the C. elegans nervous system in principle allow for the development of detailed mechanistic computational models. However, there are several challenges. We often do not have direct experimental access to important modeling details such as single-neuron dynamics and the signs and strengths of the synaptic connectivity. Further, neural activity can only be measured in a subset of neurons, often indirectly via calcium imaging, and significant trial-totrial variability has been observed. To overcome these challenges, we introduce a connectome-constrained latent variable model (CC-LVM) of the unobserved voltage dynamics of the entire C. elegans nervous system and the observed calcium signals. We used the framework of variational autoencoders to fit parameters of the mechanistic simulation constituting the generative model of the LVM to calcium imaging observations. A variational approximate posterior distribution over latent voltage traces for all neurons is efficiently inferred using an inference network, and constrained by a prior distribution given by the biophysical simulation of neural dynamics. When applied to a recent dataset, we find that connectomic constraints enable our LVM to predict the activity of neurons whose activity were withheld. We explored models with different degrees of biophysical detail, and found that models with the most realistic conductance-based synapses provide markedly better predictions than current-based synapses for this system.

heading: 1 INTRODUCTION
text: The anatomical connectivity of the entire C. elegans nervous system, including both chemical and electrical synapses, has been known for several decades [25; 23; 26]. However, well-calibrated and predictive connectome-constrained mechanistic models of this nervous system have yet to be demonstrated [10; 24; 22; 7; 5]. This is because currently available experimental data are insufficient to completely constrain computational models. First, the single-neuron and synapse dynamics are generally unknown. Second, the connectome does not directly inform the signs and strengths of individual synapses. Third, the response properties of sensory neurons are incompletely known. Further, it is unclear what level of biophysical detail is necessary to reproduce the essential computational function of the C. elegans nervous system.
We use recently available whole brain calcium imaging data to constrain the missing parameters in a connectome-constrained, biophysically detailed model of the C. elegans nervous system. We start with a simplified non-spiking passive point-neuron model of the voltage dynamics of individual neurons in the circuit. We model inputs to the neurons from electrical synapses, and model nonlinear chemical synapses with either current-based or conductance-based biophysics. The challenge of fitting such a model to data is two-fold. First, the voltage dynamics of the neurons are not directly observed, but rather indirectly measured via their slow calcium dynamics. Second, there is significant trial-to-trial variability in neural activity, suggesting strong initial state dependence in the neural responses even to the same sensory stimulus [6]. These issues can both be addressed by treating the collective voltage signals of all neurons in the nervous system as an unobserved latent variable whose dynamics are determined by the simplified connectome-constrained biophysical model with unknown neuronal and synaptic parameters.
Our connectome-constrained latent variable model (CC-LVM) of voltage dynamics of the entire C. elegans nervous system is a large-scale latent variable model with a very high-dimensional latent space consisting of voltage dynamics of 300 neurons over 5 minutes of time. Many sources say that C. elegans has 302 neurons, but two of them are not connected to the rest of the nervous system, so we only model the 300 that are [26]. The generative model for these latent variables are described by stochastic differential equations modeling the nonlinear dynamics of the network activity, and a novel differential equation describing how the calcium signals signals are generated from the voltage dynamics of individual non-spiking neurons. We developed a variational autoencoder based framework for inferring the unobserved voltage dynamics from the observed calcium dynamics of only a subset of the neurons in the nervous system. An inference network enables efficient variational inference of the latent variables, and fitting the unknown neuronal and synaptic parameters of our connectome and biophysics based generative model.
The compact nervous system of C. elegans makes it an ideal platform for systems neuroscience. The C. elegans nervous system consists of 300 neurons, divided into 118 distinct classes (often bilaterally symmetric pairs) [21; 11; 20]. These neurons can largely be categorized as sensory neurons, interneurons, and motor neurons. The majority of these neurons, about 200, are concentrated in the head, forming the brain of the worm. The synaptic connectivity of C. elegans has been mapped with electron microscopy, providing researchers with a complete connectome [25; 26].
We apply the CC-LVM to whole-brain calcium imaging data which captures 170 of the 300 neurons in multiple worms as they respond to chemosensory stimuli [27]. In principle, an accurate model of the nervous system constrained by incomplete activity measurements but complete connectivity measurements can enable accurate predictions of activity in neurons which were not recorded. We tested this hypothesis by using the CC-LVM to predict the activity of neurons which were measured, but whose activity was withheld during model training. We also used the model to predict the activity of entire worms, by holding out single trials during training. The CC-LVM predicted activity of withheld neurons and withheld worms significantly better than models unconstrained by the connectome, demonstrating the utility of the connectome even when little is known about the signs and strengths of individual connections. Further, we found that models with conductance-based synapses provide superior predictions to models with current-based synapses. CC-LVMs thus provide a new tool for connectome and activity constrained modeling of neural circuits and for discovering the appropriate level of detail of biophysical model for a given system.

heading: 1.1 PRIOR WORK
text: Previous work has been focused on creating network simulations based on anatomical connectome data, with unknown single neuron synaptic biophysics [13; 5; 24; 15; 10]. These network models provided a holistic view of the entire nervous system and were validated by comparing simulated locomotion with movements of live animals [8]. However, they were not fit against prerecorded calcium fluorescence data, so their simulated neuronal activities are not quantitatively confirmed. Recent advances in machine learning have enabled increasingly sophisticated latent variable models to uncover structure from data generated by advanced neural interfacing technologies [17]. This framework has been used in other contexts to infer neuronal voltage dynamics from high-dimensional calcium fluorescence recordings [16; 1; 19]. Despite their utility, these models do not account for connectomic data which may greatly improve their predictive power. Our work employs both LVMs and connectomic constraints to create a model that is informed by both neural population dynamics and network structure.
Recently, Bayesian models have been widely applied in neuronal datasets to infer spiking activity, neural dynamics, and connectivity [24; 1; 19]. In particular, Warrington et al. [24]. use sequential Monte Carlo (SMC) to impute the intracellular voltage potentials of C. elegans neurons from 49 recorded calcium traces [9]. Their work combined neuronal, body and calcium observation simulators in order to model locomotion [3; 14]. Their simulator produces a series of exemplar neuronal voltage traces but they do not evaluate the imputed traces of unmeasured neurons against real data so the accuracy of their inferences is unknown. Furthermore, they explore how different methods for approximate Bayesian inference affect parameter estimation on simulated data and show that their method is best. Our model uses a variational auto-encoder (VAE) to perform inference instead of SMC which is more computationally efficient because it allows direct sampling from the approximate posterior. Additionally, the CC-LVM is trained on significantly more neurons, and we validate our
model by evaluating its ability to predict the activity of neurons held out from the training data. Finally, rather than study how various inference methods can affect parameter estimation in a test autoregressive model, we search a space of generative models to optimize voltage predictions by comparing them to real calcium data.

heading: 2 CONNECTOME-CONSTRAINED LATENT VARIABLE MODEL
text: We constructed a connectome-constrained latent variable model (CC-LVM) of the C. elegans nervous system, where each node in our network represents a specific neuron in the animal. The activity of each neuron was modeled with a latent variable analogous to voltage. The dynamics of each neuron in the network is represented with a stochastic non-spiking passive leaky integrator equation with learned time constants and resting membrane potentials [4]. The neurons were coupled by both chemical and electrical synapses with learned weights [4]. We allowed these weights to be non-zero only where the connectome indicates the existence of synapses.
Given a set of learned parameters (weights, time constants, and resting membrane potentials), the dynamics of the network defines a prior distribution over neural activity trajectories. Importantly, the stochastic nature of the dynamics allows for deviations from perfectly deterministic dynamics. This allows the model to accommodate the observed variability in single-neuron dynamics. This variability has several potential sources, including the unmeasured initial states of the neurons and our incomplete knowledge of the sensory inputs driving the nervous system.
A latent variable model of this scale with a nonlinear generative model defined by the stochastic dynamics is difficult to fit. To address this challenge, we used the probabilistic inference framework of variational auto-encoders (VAE) [18] to train a black-box voltage inference network to predict a posterior distribution over neural activity trajectories. We then used this inference network to train the parameters of the CC-LVM. The resulting LVM has a biologically realistic generative model of the nonlinear neural dynamics of the C. elegans nervous system, and a black-box temporal convolutional inference network which, given sensory stimulus and calcium imaging measurements, predicts a factorized Gaussian distribution over the voltages of all the neurons in the network. Several variants of the LVM were developed: we tested models in which the synaptic connections were modeled as either current-based or conductance-based [4], and evaluated different levels of connectome constraint. The LVM was optimized with the ELBO (evidence lower bound) objective.

heading: 2.1 NETWORK MODEL WITH PASSIVE POINT-NEURON VOLTAGE DYNAMICS
text: Neurons in the C. elegans nervous system are largely non-spiking [2], so we model the voltage dynamics for these neurons as a passive point neurons with a single electrical compartment. Let v ∈ RN denote the voltages of the N neurons. The voltage vi(t) for each post-synaptic neuron i at
the time t was calculated using a first-order leaky integrator equation given by
τiv̇i(t) + vi(t) = s c i (t) + s e i (t) + v rest i + hi(t), (1)
where τi is the voltage time constant, hi is the chemosensory input provided to only the sensory neurons, sci is the chemical synaptic input, s e i is the electrical synaptic input, v rest i is the resting neuron voltage.
We studied two variations of the model, the current-based model and the conductance-based model, which differ in their formulations of chemical synaptic input sci . Since neurons in the C. elegans nervous system are largely non-spiking, we model the chemical synapses as having graded release of neurotransmitter, rather than the all-or-none quantal release seen in spiking neurons. In both models, we model the amount of neurotransmitter released as proportional to the pre-synaptic voltage vj , following softplus activation g(·) which sets a minimum voltage below which there is no synaptic release, leading to W cjig(vj(t)). We use g(·), to denote a softplus function for the rest of the paper. In our current-based model, synaptic input sci to a post-synaptic neuron i is directly proportional pre-synaptic neurotransmitter concentration:
sci (t) = N∑ j W cjig(vj(t)), (2)
where W cji represents the chemical synaptic weight between pre-synaptic neuron j and post-synaptic neuron i. W cji can be positive or negative depending on if the synaptic connection is excitatory or inhibitory. If neurons j and i are not connected, W cji is set to zero. In the conductance-based model, we model the synaptic current entering the post-synaptic neuron with more biophysical detail as
sci (t) = N∑ j (Eji − vi(t))W cjig(vj(t)). (3)
Here, the pre-synaptic neurotransmitter concentration W cjig(vj(t)) is more accurately modeled as proportional to the conductance at the post-synaptic terminal. The post-synaptic current is then given by the product of the synaptic conductance and the difference between the post-synaptic voltage vi(t) and the synaptic reversal potential Eji.
In contrast to the current-based synapse whose input is independent of the post-synaptic voltage vi(t), the conductance-based synapse model more accurately also has a dependence on the post-synaptic voltage. In addition, this model decouples the sign of the synapse from the strength. The reversal potential of a synapses Eji dictates whether a synapse is excitatory or inhibitory. A large and positive Eji corresponds to an excitatory synapse causing depolarization of the post-synaptic neuron. And an inhibitory synapse will have a negative Eji causing hyperpolarization. In this model, we can now independently train the sign of a synapse and its non-negative strength W cji, which is not easily possible with current-based synapses. In both the current-based and conductance-based models, the following equation was used to represent electrical synaptic inputs:
sei (t) = N∑ j W eji(vj(t)− vi(t)), (4)
where W eji is restricted to be non-negative and vj − vi is the potential difference between presynaptic and postsynaptic neurons. We restrictW eji = W e ij because the potential differences between electrical synapses are symmetric. To directly compare the outputs of the LVM to neural activity measurements, our model must generate calcium signals from the voltage traces. We model the calcium concentration [Ca]i of each neuron i as a first-order leaky integrator, driven by voltage-gated calcium channels with the same nonlinear current-voltage (I-V) function g(vi):
τ[Ca][Ċa]i(t) + [Ca]i(t) = g(vi(t)), (5)
where g(·) represents softplus activation, τ[Ca] is a time constant shared across all neurons. We then map calcium concentration [Ca] into the measured calcium fluorescence signals f via an affine transform with scalar αf and bias βf ,
fi(t) = α f [Ca]i(t) + β f + σf fi (t), (6)
with measurement noise represented by a noise amplitude σf and a noise term f i (t) ∼ N (0, 1).

heading: 2.2 VARIATIONAL INFERENCE
text: Note that voltage dynamics in equation 1 is a deterministic equation. We provide further flexibility to the model in order to model the significant trial-to-trial variability we observed in Figure 1c. Such variability could be caused by different initial state of neurons, or by unknown, unmeasured inputs from environment. In order to address this, we allow for the voltage dynamics to deviate from the deterministic dynamics, and introduce a Gaussian prior over these deviations. This results in the following stochastic version of the voltage dynamics,
τiv̇i(t) + vi(t) = s c i (t) + s e i (t) + v rest i + hi(t) + σ v i v i (t), (7)
Where σvi is the noise amplitude, and the standard normal noise term is v i ∼ N (0, 1). With this approach, we re-formulated the neuronal voltage as a latent variable v with the generative model given by 7 and must perform Bayesian inference to estimate their value. Exact Bayesian inference is intractable for this high-dimensional nonlinear stochastic dynamical system. For this reason, we developed a computationally efficient strategy for inferring an approximate posterior distribution P (v|f) over the latent variable v, given the measured fluorescence data f using variational Bayes in the framework of the variational autoencoder [12], see derivation in Appendix A.1.

heading: 2.3 INFERENCE NETWORK
text: The inference network Qφ(v|f ,h) generates the approximate posterior distribution of the latent variable v given the measured fluorescence data f and the sensory inputs h. It consists of multiple 1D convolutional layers that transform the input data into a Gaussian distributions vi(t) ∼ N (µi(t), σi(t)) for each neuron i = 1, ..., N at every time point t = 0, ..., D. The approximate posterior distribution Qφ(v|f ,h) is factorized as
Qφ(v|f ,h) = ∏
t=0,...,D
Qφ (v(t)|f(0), ..,f(D),h(0), ..,h(D)) , (8)
D is the length of the time window of data fed into the inference network, and φ is its parameters.

heading: 2.4 GENERATIVE MODEL
text: The generative model outputs the calcium fluorescence trace f and neuronal voltage v given sensory inputs h. It can be formulated as
Pθ(f ,v|h) = Pθ(f |v,h)Pθ(v|h), (9) where θ contains parameters of the generative model. Pθ(v|h) is a biophysically realistic connectomeconstrained network with passive point-neuron dynamics. It outputs a prior distribution over voltage v given sensory input h. We use a nonlinear sensory mapping H to transform our chemosensory input to neuronal stimulus to the sensory neurons. The stimuli is represented as a set of binary vectors hk(t), which are 1 given the presence of stimulus k at time t and 0 otherwise. Pθ(f |v,h) is another realistic model that maps the voltage v into reconstructed fluorescence f , given sensory inputs h. The voltages v are sampled from the approximate posterior distributions generated by the inference network. Details of both distributions are described in Appendix A.2.

heading: 2.5 OBJECTIVE FUNCTION FOR GENERATIVE MODEL AND INFERENCE NETWORK
text: The variational objective of our latent variable model uses ELBO, which is defined as
L = EQ(v|f ,h)[logPθ(f ,v|h)− logQφ(v|f ,h)]. (10)
The reconstructed calcium fluorescence f is then fit to the measured fluorescence data by maximizing the ELBO for individual data points. The objective function is a combination of reconstruction loss between the reconstructed and measured fluorescence traces and the KL divergence between voltage prior and posterior. (details in Appendix A.4, A.5):
L = −Ev∼Q(v|f ,h) [log (Pθ (f |v,h))] +DKL (Qφ(v|f ,h)‖Pθ(v|h)) , (11)
Optimization of the KL divergence with the reconstruction loss allows our inference network to learn the time dynamics of our neuronal voltages because the difference between the posterior and prior distributions is made small.

heading: 2.6 CONNECTOME CONSTRAINT FOR MODEL WEIGHTS
text: To build the model using chemical and electrical synapse count, Cc and Ce, as a connectome constraint, we factorized the chemical synaptic weight W cji as W c ji = α cT cjiM c ji, where α
c represents a global scaling factor which scales synapse counts to per-synapse currents or conductances depending on the model class. T cji is a binary value which indicates the connectivity, and M c ji reflects the magnitude of connection strength. The electrical synaptic weights matrix W eji is similarly defined as W eji = α eT ejiM e ji. We compared the performance of the connectome constrained models to unconstrained ones (Table 1, Appendix A.6). The anatomical connectome can be represented mathematically as two N ×N matrices: one containing the chemical synapse counts (Cc), and one containing the sizes of electrical synapses (Ce). We designed the connectome sparsity constraint by fixing the sparsity of our model to that of the connectome. For chemical connections, the model sparsity T cji = 1[C c ji > 0], and likewise for electrical connections, T e ji = 1[C e ji > 0]. For both chemical and electrical connections, the magnitudes M cji and M e ji are trainable.
We applied the connectome sparsity constraint to both current and conductance-based versions of the LVM. In the conductance-based model, we were able1 to apply the even stronger connectome count constraint: directly assigning the synapse counts to the chemical connection magnitudes (M cji = C c ji) and the synapse sizes to the electrical connection magnitudes (Meji = C e ji). Under this constraint, both the sparsity and magnitude are fixed, and only the global scaling factors αc and αe are trainable. Additionally, we evaluated another variant of this constraint, connectome count2. This constraint uses parameters of connectome count1 to initialize M cji. M c ji is then allowed to deviate from C c ji during
1Since in the current-based model, the range of W cji is (−∞,∞), and because the signs of the synaptic connections are unknown, we could not directly apply the connectome count constraint to it.
re-training. We compared our connectome-constrained models to two unconstrained models. The learned dense model used fully connected sparsity matrices T c and T e, with the magnitudes M c and M e as trainable parameters. The learned sparse model uses L1 regularization on both magnitude matrices: ‖M e‖1 and ‖M c‖1. The learned total count model preserves the total synapse count for both magnitude matrices: ‖‖M c‖1 − ‖Cc‖1‖22 and ‖‖M e‖1 − ‖Ce‖1‖22. αc and αe are 0.002 and 0.065 respectively, which is the same in the conne.count1 model. Thus, it regularizes the sum of elements in M c and M e equal to that of Cc and Ce.

heading: 3 EXPERIMENTS

heading: 3.1 DATASET
text: We applied the CC-LVM to a calcium imaging dataset in which immobilized, pan-neuronally labeled C.elegans were presented with a panel of chemosensory stimuli (2-butanone, 2,3-pentanedione, and NaCl) [27]. The activity of 170 neurons in the head of the animal were measured across 21 individuals. Each recording lasts for 5 min, with an acquisition rate of 4 Hz. In this dataset, every neuron in the worm brain could be identified in vivo using a deterministic multicolor landmark strain. Thus, it was the first experiment in which whole brain dynamics could be mapped onto the anatomical connectome at the resolution of single neurons. The connectome constraints we applied utilized the anatomical connectivity data from [25]. The complete hermaphrodite C.elegans was reconstructed using electron microscopy, and the connections between its 300 neurons were mapped. The connectome contains 3464 chemical synapses and 1031 electrical synapses. More details can be found in appendix A.6.

heading: 3.2 NEURON HOLDOUT EVALUATION
text: We hypothesized that the CC-LVM which incorporates the complete connectome but has access only to the measured activity of a subset of the neurons may nevertheless be able to accurately predict the activity of unmeasured neurons. This would constitute experimentally testable predictions from our mechanistic model. We tested this hypothesis by performing neuron holdout evaluations, withholding a single bilateral pair of measured neurons from the model during both training and testing. Note that our latent variable model always infers the voltages of all neurons regardless of how many neurons are observed in any experiment, since the latent states are the voltages of the entire nervous system. The goal of this evaluation is to test the model’s ability to predict neural activity of neurons which were never observed at any time. We fit a version of the model with a pair of neurons withheld, and then compare the model predictions for the same pair of neurons with the ground truth calcium fluorescence measurements. We repeated this evaluation for all 107 measured neuron pairs. Because many of the neurons in C. elegans are bilaterally symmetric, the symmetric pairs tend to have high correlations with each other. In order to prevent the model from predicting one neuron solely based on its sibling, we removed the symmetric neurons in groups of two. For every neuron pair, we performed the neuron holdout evaluation for both current and conductance models with different levels of connectome constraints. For these experiments, we calculated the correlation coefficient (details in Appendix A.4) between the predicted and measured calcium fluorescence magnitudes to quantify the predictive performance (Figure 3a and Table 2, more predicted traces in Appendix A.9). For each evaluation, we trained 4 models with random initialization, and reported both the mean value and standard error (SE) for each model. We found connectome sparsity and connectome count constrained models were better at predicting traces of held out neurons than unconstrained models.

heading: 3.3 WORM HOLDOUT EVALUATION
text: Another evaluation we performed was to hold out the data from a handful of individual worms, train the model on the remaining worms, and predict the activity of the neurons of the held out individuals. For each of the 9 model variants, we trained on 15 worms, and tested the model on 6 withheld worms. In each case, we trained 4 models with different random initialization. For the worm holdout traces, we calculated the correlation coefficient between the predicted and measured calcium fluorescence. We also reported ELBO, reconstruction loss and KL divergence (Appendix A.4) on tested worms to represent the fitting performance of each model. We reported the mean and standard error (SE) from models trained with 4 different random initializations. We found that models using the connectome count constraints performed better at predicting calcium traces from held out worms in Table 2.

heading: 3.4 COMPARING THE DIFFERENT MODEL VARIANTS
text: The gold standard for comparing generative models in machine learning is to compare the model evidence P (f) on a test dataset. However, this is intractable for many generative models, and so the evidence lower bound (ELBO) is used for model fitting and sometimes also for model comparison. However, the ELBO is not a tight lower bound, and so is not guaranteed to provide the same ordering over models as the model evidence P (f). In this paper, exploiting the interpretable nature of our mechanistic generative model, we developed a stricter criterion for comparing models based on the neuron holdout evaluation. Neuron holdout predictive performance is biologically meaningful, as it evaluates the ability of a model to make testable predictions about neurons which were not directly measured. We found that overall, the conductance-based CC-LVMs under connectome count constraint make the best neuron holdout predictions as quantified by the correlation coefficient. We found that model accuracies at neuron holdout prediction are not very tightly correlated with their test ELBO. While connectome count2 still achieves the best performance, the learned dense model achieves the second best performance, despite performing worst in the neuron holdout evaluation. This suggests that the learned dense model has potentially discovered an equally accurate but mechanistically incorrect way to model the calcium imaging data.
Stronger connectome constraint improves the neuron holdout predictions. We found that our CC-LVM with the connectome count constraint produced significantly better predictions than CCLVM with the connectome sparsity constraint, and both connectome constrained models outperformed unconstrained LVMs with the same dynamics (Table 2). This suggests that the neuron activity predictions are improved by stronger connectome constraint. Note that connectome count2 achieves slightly better performance than connectome count1, as its parameters get initialized from connectome count1, then are allowed to deviate during re-training. We also explored whether CC-LVMs achieve better performance due solely to the smaller number of trainable parameters. Comparing the connectomeconstrained models to the learned sparse LVM with sparsity regularization, we found that the learned sparse model performed worse. This indicates that the connectome constraint improves predictability due to its topology, not due to the number of trainable parameters.
Perhaps relatedly, we found superior predictions for interneurons and motor neurons, compared to sensory neurons. While the inputs to motor neurons and especially interneurons are largely contained in the connectome, the inputs to sensory neurons are not as well constrained by the connectome since they also respond to external stimuli. Instead, we trained a blackbox model to predict the tunings of sensory neurons to experimentally presented odors, but this likely does not fully characterize all the sensory inputs to the worm. This result additionally suggests the utility of the connectome, for neurons where the connectome provides information about the dominant sources of input.
Conductance-based models outperform current-based models. We found that conductance-based models achieve better performance in neuron holdout evaluations, despite needing to estimate the reversal potential E in addition to the non-negative strength of a chemical synapse W c. This could be for several reasons. The biophysical model of a conductance based synapse might be a more accurate description of reality. There are also two technical advantages to this model. First, conductance based synapses allow us to make the fullest use of the connectomic data by using the non-negative synapse
counts. Second, separating the optimization of the sign of a synapse from its magnitude might allow for an easier optimization of both. In contrast, the sign and magnitude of a current-based synapse are coupled in the same parameter, and we have observed strong local minima preventing the efficient optimization of the sign and the magnitude. Accordingly, we found that conductance-based models converge to a connectome closer to the counts given by the original data with less variance due to the removal of the neuron sign parameter (Appendix A.7). We found that a synapse initialized as excitatory is unlikely to change its sign after optimization. For the same reasons, we are unable to effectively use the non-negative synapse counts contained in the connectomic data to inform the magnitude of the current-based synapse while leaving the sign unconstrained.

heading: 4 DISCUSSION
text: In this work, we proposed a connectome-constrained latent variable model (CC-LVM) of the C.elegans nervous system. Our latent variable framework enables the estimation of the initial state of the network dynamics, and also captures modeling imperfections, which together enable us to fit models to calcium imaging data with significant trial-to-trial variability. In contrast to popular blackbox generative models used in LVMs in neural data modeling [17], we developed a biophysics based, connectome constrained, mechanistically detailed generative model for the neural dynamics. We evaluated two current-based and conductance-based models for synaptic connection. We also explored using connectomic data to constrain the model to different degrees of detail. Since CC-LVMs enable the prediction of neurons which were not used to fit the model, we were able to simulate and validate the process of making biologically relevant and experimentally testable predictions about neurons which have not yet been experimentally characterized.
We found the conductance-based model produces much better-correlated predictions in neuron holdout evaluation and worm holdout evaluation. And using the knowledge of the connectome can aid in silico predictions of individual neurons whose activity were not measured. We hope our mechanistic model can be used to predict the effects of perturbations or ablations, leading to a causal understanding of neural computation in neural circuits. The research in this paper was enabled by recent advances in whole brain calcium imaging in identified neurons [27]. While we have only analyzed the first datasets to be collected in this manner, we hope that a rich database of such recordings across a diversity of conditions will lead to achieving the dream of standard models of the C. elegans nervous system [8] which generalize across experiments.
Reproducibility: We trained each of our LVMs on 1 Quadro RTX 8000. Training procedures are described in Appendix A.3, A.4, A.5. We will release our software and datasets upon publication.

heading: A APPENDIX
text: A.1 VARIATIONAL INFERENCE
We used bayesian inference (VAE) to compute the posterior probability over the latent variables v (neuron voltages) given the data f (measured calcium fluorescence),
P (v|f) = P (f |v)P (v) P (f)
(12)
To solve this problem, we need to compute the model evidence,
P (f) = ∫ P (f |v)P (v)dv. (13)
However, computing P (f) is an intractable problem, so we use variational inference to approximate the value of P (v|f) [12]. Denoting the approximate posterior Q(v|f), we can formulate evidence lower bound (ELBO) as
logP (f) ≥ L = EQ(v|f)[logP (f ,v)− logQ(f |v)], (14)
which tightens the lower bound on logP (f) as L is maximized. This allows us to improve the inference network and find the posterior such that Q(v|f) ≈ P (v|f) since the Kullback-Leibler (KL) divergence, which measures the distance between the distributions Q and P , is minimized by maximizing ELBO.

heading: A.2 PRIOR AND POSTERIOR DISTRIBUTION AT A DISCRETE TIME STEP
text: Using euler discretization with a time step of ∆t, the first-order leaky integrator equation in equation 1 in the network model is
τi(vi(t+ ∆t)− vi(t))/∆t+ vi(t) = sci (t) + sei (t) + vresti + hi(t+ ∆t), (15)
Then the voltage vi(t+ ∆t) of neuron i at time step t+ ∆t is formulated as
vi(t+ ∆t) = ∆t/τi(s c i (t) + s e i (t) + v rest i + hi(t+ ∆t)− vi(t)) + vi(t), (16)
This equation indicates that the voltage state vi(t) at the current time step t is determined by a function over multiple components from the last time step t−∆t, including the past voltage state, the chemical and electric synaptic inputs as well as the sensory input
vi(t+ ∆t) = F(vi(t), sci (t), sei (t), hi(t+ ∆t)), (17)
For the posterior distribution Qφ(v(t)|f(t),h(t)), the inference network outputs µi(t) and σi(t) for every neuron i at each time step t. The posterior distribution for the voltage generated by the inference network at each time step is vi(t) ∼ N (µi(t), σi(t)), given a sequence of measured fluorescence data f(0), ...,f(D∆t) and sensory input h(0), ...,h(D∆t) with total number of steps D + 1. Then we used the reparameterization trick to sample the voltage vi(t) from this posterior distribution, to generate a sequence of voltage sample from vi(0), ..., vi(D∆t).
For the prior distribution, we used the voltage samples vi(t) from the posterior distribution at the last step as the past voltage state in F to calculate the prior distribution Pθ(vi(t + ∆t)|hi(t + ∆t)) at time step t+ ∆t. The prior distribution is defined as N (F(vi(t), sci (t), sei (t), hi(t+ ∆t)), σvi ). σvi is the constant standard deviation of the prior distribution that is trainable during optimization. At the initial time step 0, we define another trainable parameter µi(0) to represent the mean of the voltage prior distribution.
Similarly, we also used euler integration to discretize the first-order leaky integrator equation 5. The calcium concentration [Ca]i(t+ ∆t) of neuron i at time step t+ ∆t is given by
[Ca]i(t+ ∆t) = ∆t/τ[Ca](g(vi(t))− [Ca]i(t)) + [Ca]i(t), (18)
We use a sequence of voltage samples vi(t), t = 0, ..., D∆t as the inputs to this formula to get a sequence of calcium concentration [Ca]i(1), ..., [Ca]i(D∆t). The initial state of concentration is estimated from the measured fluorescence trace at time step 0 using [Ca]i(0) = (Fi(0)− βf )/αf . equation 6 in the generative model Pθ generates the output distributions Pθ(f(t)|v(t),h) of fluorescence trace at time step t as N (αf [Ca]i(t) + βf , σf ). The values αf , βf , σf are trainable fluorescence parameters.

heading: A.3 MODEL ARCHITECTURE
text: The calcium fluorescence data is given in a matrix of dimensions (N,T ) where N is the total number of neurons in the connectome and a matrix of dimensions (3, T ) which is 3 one-hot encoding vectors that are 1 when the chemosensory stimulus is present and 0 otherwise. Each of the 3 vectors represent one of the three stimulus presented: NaCl, 2-butanone and 2,3-pentanedione.
The calcium fluorescence matrix is passed to the input layer of our inference networkQ(v|f ,h) which is a 1D convolutional layer. Then the output of the input layer gets passed through an upsampling layer, a ReLU nonlinearity, then another 1D convolutional layer and then another upsampling layer. The first convolutional layer has 2N input channels, 2N output channels and a kernel size of 11 and the second has 2N input channels, 2N output channels and a kernel size of 21. The first upsampling layer upsamples by a factor of 4 and the second upsamples by a factor of 10 to account for the fact that we simulate the network at a rate 40x faster than the recording rate of our data. The sensory input is passed through a 1D convolutional layer with N input channels, 2N output channels and a kernel size of 21. This is passed through a ReLU and then concatenated to the output of the 2nd convolution on our calcium data. The merged matrix is then 1D convolved again with 2N inputs, 2N outputs and a kernel size of 41. This is one of the final outputs of our inference network and gives a (N,T ) matrix that represents the means of our posterior distribution. To generate the standard deviation of our factorized posterior distribution, we pass our merged data through another convolutional layer with the same dimensions but different parameters. We pass that output through a nonlinear affine transform (using a softplus nonlinearity) to rescale the output to an appropriate level of noise. The parameters of this inference network are all the weights and biases of the convolutions and upsampling layers described above, which we denote by φ in our paper.
The input to the generative model is the same sensory input that we pass to our inference network. The first layer of our generative model is a single layer nonlinear transform that converts our 3 dimensional input to an N dimensional input that is 0 to any non-sensory neurons. The next layer is our recurrent neural dynamics, given by equation 1, which passes information to itself in the next time step based on the connectivity given by the connectome. The voltage simulated by each neuron in the connectome is passed to our calcium simulation layer, equation 5. The final layer is then a calcium to fluorescence affine transformation, equation 6. The output of this final layer is then compared to the ground truth fluorescence data to calculate the reconstruction loss. Below, we have included some psuedocode describing the architecture of our network
I n f e r e n c e _ N e t w o r k ( c a l c i u m _ f l u o r , mi s s ing_da ta_mask , s e n s o r y _ i n p u t ) : c a l c i u m _ i n p u t = c o n c a t e n a t e ( c a l c i u m _ f l u o r , m i s s i n g _ d a t a _ m a s k ) conv1_ou t = r e l u ( conv1 ( c a l c i u m _ i n p u t ) ) up1_out = upsample ( conv1_ou t ) conv2_ou t = r e l u ( conv2 ( up1_out ) ) up2_out = upsample ( conv2_ou t )
s e n s o r y _ c o n v _ o u t = r e l u ( conv3 ( s e n s o r y _ i n p u t ) ) m e r g e d _ c a l c i u m _ s e n s o r y = c o n c a t e n a t e ( up2_out , s e n s o r y _ c o n v _ o u t )
m e a n _ l a t e n t _ n e u r o n _ v o l t a g e = conv4 ( m e r g e d _ c a l c _ s e n s o r y )
s t d _ l a t e n t _ n e u r o n _ v o l t a g e = s o f t p l u s ( conv5 ( m e r g e d _ c a l c i m _ s e n s o r y ) )
s a m p l e _ l a t e n t _ n e u r o n _ v o l t a g e = m e a n _ l a t e n t _ n e u r o n _ v o l t a g e
+ rand_norm * s t d _ l a t e n t _ n e u r o n _ v o l t a g e
r e t u r n s a m p l e _ l a t e n t _ n e u r o n _ v o l t a g e
Gene ra t i ve_Mode l ( s a m p l e _ l a t e n t _ n e u r o n _ v o l t a g e , s e n s o r y _ i n p u t ) :
# E q u a t i o n 1 n e u r o n _ v o l t a g e _ d y n a m i c s = l e a k y _ i n t e g r a t o r _ c o n n e c t o m e _ d y n a m i c s ( s a m p l e _ l a t e n t _ n e u r o n _ v o l t a g e , s e n s o r y _ i n p u t )
# E q u a t i o n 5 c a l c i u m _ c o n c e n t r a t i o n _ d y n a m i c s = l e a k y _ i n t e g r a t o r _ c a l c i u m _ m o d e l ( n e u r o n _ v o l t a g e _ d y n a m i c s )
# E q u a t i o n 6 f l u o r e s c e n c e _ t r a c e = n o n l i n e a r _ a f f i n e _ t r a n s f o r m ( c a l c i u m _ c o n c e n t r a t i o n _ d y n a m c i s )
r e t u r n f l u o r e s c e n c e _ t r a c e

heading: A.4 FORMULATIONS FOR OBJECTIVE FUNCTIONS AND METRICS
text: Here we provided the mathematical formulations for the objective functions and metrics we used in this work.
The evidence lower bound (ELBO) is a combination of reconstruction log likelihood and the negative KL divergence. The reconstruction log likelihood, given the predicted fluorescence trace f̂i(t) of neuron i at time step t as output distribution N (f̂i(t), σf ). And f̂i is formulated as
f̂i(t) = α f [Ca](t) + βf . (19)
Given the measured fluorescence trace fi(t), number of neurons N , and total number of time steps D + 1, we calculated the reconstruction negative log likelihood as the reconstruction loss
Lrecon = − D∆t∑ t=0 (2π(σf )2)−N/2 exp
( − 1
2(σf )2 N∑ i=1
(fi(t)− f̂i(t))2 )
(20)
For KL divergence between voltage prior and posterior distribution. We represent the posterior distribution generated from the inference network as N (µi(t), σi(t)), for neuron i at each time step t. We then sample vi(t) from the posterior distribution. According to equation 17, we represent the mean v̂i(t) of prior distribution N (v̂i(t), σvi ) as
v̂i(t) = F(vi(t), sci (t), sei (t), hi(t+ ∆t)), (21)
where σvi is a trainable parameter.
Then we use a closed form of KL divergence defined between two univariate Gaussian distributions
DKL = D∆t∑ t=0 N∑ i=1 ( log ( σvi σi(t) ) + (σi(t)) 2 + (µi(t)− v̂i(t))2 2(σvi ) 2 − 1 2 ) (22)
We minimize the objective function L = −ELBO during optimization given by
L = Lrecon +DKL (23)
In the worm holdout evaluation, for each neuron pair, we calculate those metrics including ELBO, reconstruction loss and KL divergence on each holdout worm, and in total 6 worms. We also train
each model setting under the same constraint with different initialization for 4 times. We calculate the mean and standard error of ELBO, reconstruction loss and KL divergence based on 4× 6 = 24 samples for each model as our metrics.
For both neuron and worm holdout evaluation, we calculate the correlation coefficients between the measured fi(t) and reconstructed fluorescence traces f̂i(t), the correlation coefficient for each neuron is simply
ri =
D∆t∑ t=0 ( f̂i(t)− ¯̂fi(t) )( fi(t)− f̄i(t) ) √√√√D∆t∑
t=0
( f̂i(t)− ¯̂fi(t) )2 D∆t∑ t=0 ( fi(t)− f̄i(t) )2 . (24)
In the neuron holdout evaluation, for each neuron pair, we train a model in which the ground truth of said neuron pair was removed from the training data. In total there are 107 pairs, and 170 measured neurons in total. We also trained 4 times with different initialization for each neuron pair and each of 8 connectome constraints for a total of 8× 107× 4 = 3424 different models. We calculate the mean and standard error of ri based on 4× 170 = 680 samples for each model as our metrics. In the worm holdout evaluation, we calculate this metrics on each holdout worm, and in total 6 worms. Each worm has 170 measured neurons. We also train each model setting under the same connectomic constraint with different initialization 4 times. We calculate the mean and standard error of ri based on 6× 170× 4 = 4080 samples for each model as our metrics.

heading: A.5 GRADIENT-BASED OPTIMIZATION
text: We used PyTorch to minimize the ELBO using gradient descent with respect to the parameters of both the generative and inference models,
For LVMs which used connectome count constraint with αc and αe trainable. The objective function is
L = L(τ, τ[Ca], αc, αe, αf , βf , σv, σf ,E,MLP, φ) (25)
For LVMs which used other constraints including connectome count and learned dense, learned sparse, where M c and M e trainable. The objective function is
L = L(τ, τ[Ca],M c,M e, αf , βf , σv, σf ,E,MLP, φ) (26)
where φ is all the parameters of the inference network. We used a clamp in PyTorch to constrain the weight magnitude M c and M e to be non-negative in several models. We restrict the M e to be symmetric by defining it as the sum of an upper triangular matrix and its transpose. Our simulation time step is 6.25ms. Each window size D lasts for 7.5s with 30 imaging steps. We used weights of 1 for both the reconstruction loss and KL divergence. We used Adam with a learning rate scheduler to perform the optimization. We used a initial learning rate of 3e−4, with a learning rate schedular with step size of 50, and gamma of 0.5. We also set a gradient clip value of 1. Each model was trained for 300 epochs in which each epoch is one full pass through all the training data.
In the main paper, we use the variable θ to denote the trainable parameters of our generative model. The parameters in θ are given by
θ = {τ, τ[Ca],M c,M e, αf , βf , σv, σf ,E, H} (27)
where τ and τ[Ca] are the time scales of the neuron and calcium dynamics. M c and M e are the magnitudes of neuronal connections for chemical and electrical synapses respectively. αf and βf are the scaling and bias parameters of the fluorescence affine transform. σv and σf are the noise magnitudes in the neural dynamics equation 1 and the fluorescence equation 6. The variable H represents the parameters of the sensory mapping that transforms chemosensory data into neural stimulus.
The variable φ denotes all the trainable parameters of the inference network. These are just all the weights and biases of the 5 different convolutional layers described in section A.3.

heading: A.6 CONNECTOME CONSTRAINTS FOR MODEL WEIGHTS
text: We constructed CC-LVMs on both current-based and conductance-based synapses with different levels of connectomic constraints. Given anatomical connectome data with two N ×N matrices: chemical synapse counts (Cc), and sizes of electrical synapses (Ce). The chemical synaptic weight W cji as W cji = α cT cjiM c ji, where α
c represents a global scaling factor. T cji is a binary value which indicates the connectivity, and M cji reflects the magnitude of connection strength. The electrical synaptic weights matrix W eji is similarly defined as W e ji = α eT ejiM e ji. We provide more implementation details of each constraint below,
conductance-based + connectome count1 We directly assign the synapse counts to the chemical connection magnitudes (M cji = C c ji) and the synapse sizes to the eletrical connection magnitudes (Meji = C e ji). Under this constraint, both the sparsity T c ji and T e ji are fixed to the connectome sparsity, magnitude M cji and M e ji are fixed to the connectome synapse count, and only the global scaling factors αc and αe are trainable. M cji and M e ji are restricted to be non-negative.
conductance-based + connectome count2 We initialize all model parameters with the parameters from the trained conductance-based + connectome count1 model. And then the sparsity T cji and T eji are fixed to the connectome sparsity, global scaling factors α
c and αe are fixed to 0.01. And now magnitude M cji and M e ji are trainable during the re-training. M c ji and M e ji are restricted to be non-negative.
conductance-based + connectome sparsity For chemical and electrical connections, the model sparsities are T cji = 1[C c ji > 0] and T e ji = 1[C e ji > 0] respectively. For both chemical and electrical connections, the magnitudes M cji and M e ji are trainable. The sparsity T c ji and T e ji are fixed to the connectome sparsity and global scaling factors, αc and αe, are fixed to 0.01. M cji and M e ji are restricted to be non-negative.
conductance-based + learned total count For chemical and electrical connections, the model is fully connected with the sparsity T cji = 1 and T e ji = 1 respectively. For both chemical and electrical connections, the magnitudes M cji and M e ji are trainable. The model is constrained to preserve the same total synapse count as connectome count ‖Cc‖ for both magnitude matrices adding ‖‖M c‖1− ‖Cc‖1‖22 and ‖‖M e‖1 − ‖Ce‖1‖22 in the objective function. αc and αe are defined as 0.002 and 0.065, which are the same as the trained conne.count1 model, so that regularizes the sum of elements in M c equal to that of Cc, and likewise for M e and Ce. The sparsity T cji and T e ji are fixed to 1 and M cji and M e ji are restricted to be non-negative.
conductance-based + learned sparse For chemical and electrical connections, the model is fully connected with the sparsity T cji = 1 and T e ji = 1 respectively. For both chemical and electrical connections, the magnitudes M cji and M e ji are trainable. Meanwhile, it uses an additional L1 regularization on both magnitude matrices: ‖M c‖1 and ‖M e‖1, so that the fraction of nonzero elements in M c is small, and likewise for M e. The sparsity T cji and T e ji are fixed to 1 and global scaling factors αc and αe are fixed to 0.01. M cji and M e ji are restricted to be non-negative.
conductance-based + learned dense For chemical connections and electrical connections, the model is a fully connected with sparsity T cji = 1 and T e ji = 1 respectively. For both chemical and electrical connections, the magnitudes M cji and M e ji are trainable. The sparsity matrices, T c ji and T e ji are fixed. Global scaling factors αc and αe are also fixed to 0.01. M cji and M e ji are restricted to be non-negative.
current-based + connectome sparsity For chemical and electrical connections, the model sparsity is T cji = 1[C c ji > 0] and T e ji = 1[C e ji > 0] respectively. For both chemical and electrical connections, the magnitudes M cji and M e ji are trainable. The sparsity T c ji and T e ji are fixed to the connectome sparsity and global scaling factors, αc and αe, are fixed to 0.01. Meji are restricted to be non-negative, while the sign of M cji is not restricted.
current-based + learned sparse For chemical and electrical connections, the model is fully connected with sparsity T cji = 1 and T e ji = 1. For both chemical and electrical connections, the magnitudes M cji and M e ji are trainable. For both chemical and electrical connections, the magnitudes M c ji and Meji are trainable. Meanwhile, it uses an additional L1 regularization on both magnitude matrices:
||M c||1 and ||M e||1, so that the fraction of nonzero elements in M c is small, and likewise for M e. The sparsity T cji and T e ji are fixed to 1, global scaling factors, α
c and αe, are fixed to 0.01. Meji are restricted to be non-negative, while the sign of M cji is not restricted.
current-based + learned dense For chemical and electrical connections, the model is fully connected with T cji = 1, and T e ji = 1. For both chemical and electrical connections, the magnitudes M c ji and Meji are trainable. The sparsity T c ji and T e ji are fixed to 1 and global scaling factors α
c and αe are fixed to 0.01. Meji are restricted to be non-negative, while the sign of M c ji is not restricted.

heading: A.7 MODEL WEIGHT CONSISTENCY ANALYSIS
text: We used the correlation coefficient to quantify the inter-consistency of models trained with different random initializations, as well as the connectome-consistency, which is the correlation coefficient between model weights and connectome counts.
In particular, for the connectome consistency, we compare the consistency of our model weights and connectome synapse count. We calculate the correlation coefficient ri between the absolute value of model weights |Wij | and connectome synapse count Cij along each row i. Then it is averaged across N rows and 4 models with identical constraints but different random initializations to get r̄i as our metric for connectome-consistency,
ri =
N∑ j=1 ( |Wij | − ¯|Wij | )( Cij − C̄ij ) √√√√ N∑
j=1
( |Wij | − ¯|Wij | )2 N∑ j=1 ( Cij − C̄ij )2 . (28)
For inter-consistency, we evaluate how model weights are affected by random initialization after training. We calculate the correlation coefficient ri between the weights W 1ij of model 1 and weights W 2ij of model 2 along each row i. Model 1 and model 2 have the same constraints but different random initializations. The correlation coefficient between them is
ri =
N∑ j=1 ( W 1ij − W̄ 1ij )( W 2ij − W̄ 2ij ) √√√√ N∑
j=1
( W 1ij − W̄ 1ij )2 N∑ j=1 ( W 2ij − W̄ 2ij )2 . (29)
This coefficient, ri, is then averaged across N rows and each combination of 2 models within the 4 models with random initialization to get r̄i as our metric for inter-consistency.
As shown in Table 3, the conductance-based model parameters are more consistent under different random initializations. Optimization of the conductance-based model is easier because the chemical
weights W c of the conductance-based synapses are always positive and we do not need a trainable sign variable. Meanwhile, the model parameters are better constrained with the connectome sparsity or counts, which enables it to be consistent across different initialization.

heading: A.8 PREDICTIVE PERFORMANCE FOR WORM HOLDOUT V.S. NEURON HOLDOUT
text: We firstly demonstrate the violin plot for worm holdout evaluation in Figure 4. We find that the results from worm holdout is different from those from neuron holdout in Figure 3. Here we noticed that all 9 LVM variants have very close performance in different categories. And our connectome-constrained models do not outperform other unconstrained baselines.
Furthermore, we also study how the predictive performance of worm holdout and neuron holdout are correlated. In Figure 5, we show the scatter plot of predictive performance for worm holdout v.s. neuron holdout. Each point represents one neuron from the same worm. And we also report the correlation coefficient between the predictive performance for worm holdout and neuron holdout for each LVM variant. Interestingly, we find our model using connectome constraint achieves highest correlation, which indicates highest consistency for neuron sets which gets better performance in both evaluations.

heading: A.9 PREDICTED HOLDOUT TRACES WITH VARIOUS CONSTRAINTS
text: For each neuron, we performed both neuron holdout evaluation and worm holdout evaluation for current-based and conductance-based models with different levels of connectomic constraint.
For neuron holdout evaluation, we plot additional predicted traces from 15 representative neurons, including 12 with good performance (top), and 3 with worse performance (bottom) in Figure 6, model with connectome count constraint outperforms other baselines.
For worm holdout evaluation, we show the results on the same 15 neurons from 9 LVMs in Figure 7, all models achieves comparable results.

heading: A.10 PREDICTED TRACES AND VOLTAGE LATENT SPACE FOR NEURON HOLDOUT EVALUATION
text: For neuron holdout evaluation, we also show the predicted traces and voltage latent space for all measured neurons. We show 189 neurons in Figure 8, Figure 9 and Figure 10. In each subplot, we reported the measured fluorescence trace ground truth, predicted fluorescence trace ground truth, and the mean and standard deviation (std) of voltage latent space (posterior distribution) from the inference network.

heading: A.11 PREDICTED TRACES AND VOLTAGE LATENT SPACE FOR WORM HOLDOUT EVALUATION
text: For worm holdout evaluation, similarly, we also show the predicted traces and voltage latent space for all neurons. We show 300 neurons from a train worm in Figure 11–Figure 15. Meanwhile, we also show 300 neurons from a test worm in Figure 16–Figure 20. In each subplot, we reported the measured fluorescence trace ground truth (for recorded neurons), predicted fluorescence trace ground truth, and the mean and standard deviation (std) of voltage latent space (posterior distribution) from the inference network.
For the train worm, we find most of the neurons have good predictive performance compared to the measured ground truth. That indicates optimization is saturated on train set. However, test worm has a relatively worse predictive performance, that indicates the generalization gaps between each individual trial (animal). While some neurons still achieve good results, that indicates our model is still generalized to unseen individuals to some degree.
Interestingly, for both train and test worms in the worm holdout evaluation, we find those unmeasured neurons still have reasonable patterns. That indicates our model successfully avoids the degeneracy of solutions for those unmeasured neurons.

Review summaries:
The authors propose a biologically constrained latent linear dynamical model of the C. elegans nervous system. They use connectomic information (including chemical vs electrical synapses) to constrain connections between units during inference. They fit the model to calcium imaging data from whole C. elegans using a variational approach. They find that biological constraints improve model performance and validity using both withheld neuron and across worm metrics.

The paper develops a latent variable model with biologically meaningful parameters for the worm *C. elegans* whole brain neural traces. Given the connectivity between the neurons, the neural traces are modeled as stochastic leaky integrators with either conductance or current based synaptic input model. Given the voltage traces, the observed calcium signals are modeled as first order leaky integrators added by noise to account for the observation noise. External input is nonlinearly transformed and fed into the neurons to account for the stimulus-dependent modulation of neural activity. Multiple connectivity models (fully trainable, fully trainable with $L_2$ regularization, trainable weight magnitudes, trainable global scale) are used to investigate whether biological connectome constraints help with neural activity predictions or not. Variational inference is used to infer a posterior distribution over neural activity trajectories with factorized Gaussian distribution as the variational family. Results on neuron-hold-out and worm-hold-out experiments suggest that connectome constraints and conductance-based modeling improves model prediction over unseen neurons and worms.

The current manuscript presented a connectome constrained latent variable model for whole-brain calcium imaging of C.elegans nervous system. The whole-brain calcium imaging of C.elegans was collected while C.elegans was undergoing chemosensory testing, and the dataset has been published recently. In the current study, the authors aimed to present a model that could predict the single-neuron and the single-trial activity of this dataset. Specifically, the activity of each neuron of the whole brain imaging was modeled by latent variable analogous to the voltage of one unit in the model network. The connection between model network units was constrained by the connectome. The authors showed that the connectome constraints significantly improved the prediction power of the model, in predicting the activity of missing units, as well as predicting the single-trial activity of hold-out warms. Overall, the authors have provided a clear description of the model and carried out a good amount of experiments to evaluate different model variants.