Developing a 3-D computational model of neurons in the central amygdala to understand pharmacological targets for pain

Neuropathic and nociplastic pain are major causes of pain and involve brain areas such as the central nucleus of the amygdala (CeA). Within the CeA, neurons expressing protein kinase c-delta (PKCδ) or somatostatin (SST) have opposing roles in pain-like modulation. In this manuscript, we describe our progress towards developing a 3-D computational model of PKCδ and SST neurons in the CeA and the use of this model to explore the pharmacological targeting of these two neural populations in modulating nociception. Our 3-D model expands upon our existing 2-D computational framework by including a realistic 3-D spatial representation of the CeA and its subnuclei and a network of directed links that preserves morphological properties of PKCδ and SST neurons. The model consists of 13,000 neurons with cell-type specific properties and behaviors estimated from laboratory data. During each model time step, neuron firing rates are updated based on an external stimulus, inhibitory signals are transmitted between neurons via the network, and a measure of nociceptive output from the CeA is calculated as the difference in firing rates of pro-nociceptive PKCδ neurons and anti-nociceptive SST neurons. Model simulations were conducted to explore differences in output for three different spatial distributions of PKCδ and SST neurons. Our results show that the localization of these neuron populations within CeA subnuclei is a key parameter in identifying spatial and cell-type pharmacological targets for pain.

Each of the 13,000 neurons has ten variables ( Table A1) defining its properties and behaviors. During initialization, each neuron is assigned an ( , , ) location within the CeA spatial domain and a protein-expression type (Type) equal to either PKC or SST. A neuron's location and protein-expression type remain constant throughout a simulation. Each neuron has a firing frequency pattern (Freq) equal to regular spiking (RS), late firing (LF), or spontaneous (Spont). Each neuron has three variables ( , ! , " ) related to "damage". The damage variable (d) represents the percentage of total damage accumulated by a neuron in response to external stimulation. The rate at which a neuron accumulates damage during stimulation depends on its damage latency period ( ! ) and sensitivity ( " ). Additionally, each neuron has three variables related to its connectivity with other neurons within the CeA. Two of these variables are the neuron's number of incoming connections (Num-In-Link) and number outgoing connections (Num-Out-Link). The third variable is the neuron's inhibition status (Inhibited), which is a binary variable indicating whether the neuron is inhibited or not. Lastly, each neuron has a firing rate (Fr) describing the frequency in hertz (spikes per second) of the neuron's action potentials.
Each connection (i.e., directed link) in the neural network has three variables (Table A2) describing how inhibitory signals are transmitted between its endpoints. Each connection has a transmitting endpoint ( 1) equal to the ID of the agent that is the source of the signal and a receiving endpoint ( 2) equal to the ID of the agent that is the destination of the signal. Each connection has a signal strength (str) equal to the firing rate of the neuron associated with its transmitting endpoint. When used, the neural network is established during the model's initialization and does not change during a simulation. A user can choose not to use the neural network by switching the 'Network' toggle on the user interface to "Off".

Global Variables and Input Data
Global variables #$% and #$% control the size of the neural network.
#$% specifies the maximum number of outgoing connections assigned to each neuron.
#$% specifies the search radius used by a neuron when forming links with other neurons.
Both #$% and #$% are non-negative values that depend on the neuron's type (PKC or SST) and can be specified on the model's interface.
The timing, duration, and magnitude of stimulation (measured in pA) must be input as a file consisting of integer values, ranging from 0 to 220. During initialization, the values in the file are read one at a time and stored as a global vector, , where the & represents the stimulation (pA) applied during the th time step. In our simulations, the current ranges from 120 pA to 220 pA.
Lastly, the global variable '(' tracks the cumulative number of time steps during which stimulation is greater than or equal to 120 pA. During initialization, this variable is set to 0, and is increased by 1 during times steps in which & ≥ 120.

Process overview and scheduling
The following processes occur in the order listed during each simulation of the ABM.
1. Input data is provided by user.
2. Initialization of model: i. Read in patch coordinates associated with CeC, CeL, and CeM from a file and create spatial domain representing the CeA.
ii. Create 13,000 agents representing neurons in the CeA.
iii. Assign attributes to all neurons.
iv. Create network of directed links between neurons (if network is turned on).
v. Create vector specifying stimulation history. v. Use network to send inhibitory signals between neurons. If a neuron is inhibited, its firing rate is set to zero.
vi. Update all system-level observations.

Emergence
The primary emergent feature of the model is a measure of nociceptive output from the CeA that evolves over time and in response to noxious stimulation. In the model, nociceptive output is measured as the difference between the cumulative firing rates of pro-nociceptive PKC neurons and the cumulative firing rates of anti-nociceptive SST neurons in the CeA.
During each model time step, individual neuron properties are updated and then a system-level measure of nociceptive output is calculated.
Other emergent features of the model include the size of the neural network (i.e. total number of links), average number of incoming links for PKC and SST neurons, and average number of outgoing links for PKC and SST neurons.

Adaptation
A damage accumulation model is used in the ABM to allow neurons to adapt to a sensitized state over time and during stimulations measuring 120 pA or higher. Neurons start with zero damage and accrue damage at individual rates and only during time steps in which & ≥ 120. As a neuron's damage increases, the neuron transitions towards a sensitized state and its firing rate is adjusted accordingly. A neuron is considered fully sensitized when its damage variable has reached maximum value ( = 100).
Based on lab data, the firing frequency of some SST neurons changes with injury.

Interaction
Interaction occurs in the model when a neuron is inhibited by one or more other neurons.
If the sum of the inhibitory signals transmitted to a neuron during a single time step exceeds the inhibition threshold, the neuron is inhibited (i.e., & = 0) during that time step. The inhibition threshold is set to 15 Hz in all of our simulations; however, users can modify the inhibition threshold using a slider on the model interface.

Stochasticity
Stochasticity occurs at multiple time points during the model's initialization and procedures.
During initialization, each neuron is assigned a random location within the CeA or within a specific subnucleus (CeC, CeL, or CeM) if 'Non-Uniform Distribution' is selected (Section 6).
Values of damage parameters tS and tL are randomly determined for each neuron using a uniform probability distribution with ranges displayed in Table A1. Neural connectivity is achieved by creating a network of directed links between neurons using a stochastic algorithm (Section 7.1).
During each time step, a neuron's firing rate is stochastically updated using a weighted sum of values selected from truncated normal distributions (Section 7.2). Due to the stochastic nature of the model, each model simulation in our study is repeated 100 times and the mean, standard deviation, and confidence interval for all emergent output is reported.

Initialization
The first step in the model's initialization is the creation of CeA spatial domain. The CeA has three subnuclei: CeC, CeL, and CeM. Premade files containing the 3-dimensional coordinates of patches associated with each subnucleus are read into the model. Patches within each subnucleus are labeled and assigned the same color (CeC = red, CeL = blue, CeM = green). Next, 13,000 neurons with the variables in Table A1 are created and assigned a location within the CeA according to either a uniform or non-uniform distribution. If the user selects 'Uniform Distribution' on the model's interface, each neuron will be randomly assigned a location within the CeA. If the user selects 'Non-Uniform Distribution' on the model's interface, then they must also specify the percentage of PKC and SST neurons within each subnucleus of the CeA using sliders on the interface. The appropriate number of PKC and SST neurons will then be randomly assigned to a location each subnucleus. Once all neurons are created and assigned a location, a network of directed links is created between the neurons.

Creation of Network
The following procedures are applied to each neuron, one at a time, to create a network of unidirectional links between neurons. The selected neuron is identified as the transmitting neuron.
The transmitting neuron's number of outgoing links (Num-Out-Link) is set to zero. An empty list is created to track the agent IDs of all neurons receiving a link from the transmitting neuron.
While the transmitting neuron's Num-Out-Link is less than the maximum number of links allowable ( #$% ), the neuron will attempt to make new links with other nearby neurons.
To do this, the model randomly selects a receiving neuron that is not already connected to the transmitting neuron (i.e., agent ID of receiving neuron is not the list) and is located within the maximum distance ( #$% ) of the transmitting neuron. If such a neuron exists, a directed link is created from the transmitting neuron to the receiving neuron, the connectivity variable (Num-Out-Link) of the transmitting neuron is incremented by 1, and the agent ID of the receiving neuron is added to the list. If no such neuron exists, no additional links are created from this transmitting neuron. The network algorithm terminates when all neurons have attained their maximum number of outgoing links or no more suitable connections can be formed.

Update of Damage Variable for Individual Neurons
where & is the value of the damage variable at time step , " is the length of the neuron's sensitization period, and ! is the length of the neuron's latency period.

Update of firing rates for individual neurons
During each time step, the firing rates of all late firing and regular spiking neurons are stochastically updated using the equation where & is the neuron's damage level at time step , and and are random variables  Table A3. Equation (2) is linear combination of and such that when the neuron has no damage ( = 0), the firing rate of the unsensitized neuron is updated using the variable only. When damage reaches its maximum value ( = 100), the firing rate of the sensitized neuron is updated using the variable only.

Application of Network to Inhibit Neurons
After the firing rates of all PKC and SST neurons are updated, the neural network is used to transmit inhibitory signals between neurons in the ABM. The strength of an inhibitory signal transmitted through a directed link is equal the firing rate of the neuron on the transmitting end.
All PKC and SST neurons are evaluated one at a time and in a random order. For each neuron, if the total strength of all incoming signals is greater than or equal to 15 Hz, the neuron is inhibited (i.e., firing rate set to zero). If the total strength is less than 15 Hz, the neuron's firing rate does not change.

Calculation of Nociceptive Output
At the end of each time step (3) where & is a neuron's damage and & is a neuron's firing rate during time step . The first summation in equation (3) represents the weighted sum of firing rates over all PKC neurons that are either LF or RS. Each firing rate is weighted by the damage level of the corresponding neuron. As such, PKC neurons do not contribute to nociceptive output when damage is zero (i.e. pre-injury), but gradually contribute to output as sensitization occurs. When all PKC neurons have become sensitized (i.e., & = 100), each LF and RS PKC neuron contributes its firing rate to the model's output. The second summation in equation (3) represents the sum of firing rates over all SST neurons that are either late firing or regular spiking. It is assumed that SST neurons contribute to nociception at all time steps regardless of damage. Due to their antinociceptive properties, SST neurons are assumed to have a negative impact on the model's output.

Implementation
The model was coded in NetLogo3D Version 6.2.0 (Wilensky 1999). This software has a unique programming language and customizable interface that is designed specifically for ABM development and implementation. We designed a GUI for our ABM that allows a user to easily modify parameters values, network settings, and the stimulation history. The Netlogo3D code and input files for simulating the ABM can be found in an Open Science Framework public repository (https://osf.io/xnrqa/, doi 10.17605/OSF.IO/XNRQA). For the results presented in this manuscript, we used BehaviorSpace within NetLogo3D to automate batches of 100 replicate simulations for each scenario. All graphical and statistical analyses of model output were conducted in R (R Core Team 2016).   Table A3: Parameters defining the probability distributions for random variables X (unsensitized firing rate) and Y (sensitized firing rate) in equation (2). Both X and Y have truncated normal distributions with mean , standard deviation , minimum value , and maximum value .