Ping tutorial updates (#476)

* Add citation links

* Add alternative method for creating edges from adjacency matrix

* New PING network illustration

Thanks to Annabel for creating a much cleaner version!

* Final tweaks to address direct comments

* Updating documentation for PING blocks

Missing docstrings for API creation and Scott's right we should have that for tutorials

* Adding bullet points to the intro

* Suggestions are pretty much covered now

Still missing API links though - will discuss further before merging.

* Just seeing if this compiles the docs

* and does this compile

* try adding a note

* add offset to last line

---------

Co-authored-by: haris organtzidis <[email protected]>

docs/src/tutorials/ping_network.jl

@@ -1,8 +1,14 @@
 # # Pyramidal-Interneuron Gamma network
 # # Introduction
-# This tutorial provides a simple example of how to use the Neuroblox package to simulate a pyramidal-interneuron gamma (PING) network. 
+# This tutorial provides a simple example of how to use the Neuroblox package to simulate a [pyramidal-interneuron gamma (PING) network](https://direct.mit.edu/neco/article-abstract/17/3/557/6926/Effects-of-Noisy-Drive-on-Rhythms-in-Networks-of?redirectedFrom=fulltext). 
 # These networks are generally useful in modeling cortical oscillations and are used in a variety of contexts.
 # This particular example is based on Börgers, Epstein, and Kopell [1] and is a simple example of how to replicate their initial network in Neuroblox.
+# This tutorial will cover the following concepts:
+# - Creating populations of neurons for rapid simulation
+# - Creating input current sources drawn from user-specified distributions
+# - Creating a directed graph of connections between neurons
+# - Compiling the graph into a system of ODEs and simulating the network
+# - Visualizing the results with raster plots
 
 # ## Conceptual definition
 # The PING network is a simple model of a cortical network that consists of two populations of neurons: excitatory and inhibitory.
@@ -54,6 +60,7 @@ g_EI = 0.8; ## excitatory-inhibitory connection weight, manually tuned from valu
 # Finally, setup the driving currents. All neurons receive a base external current, and the inhibitory and driven excitatory populations receive a second external stimulus current.
 # The undriven excitatory neurons receive a small addition to the base current in lieu of the stochastic current in the original implementation.
 # There is also an external inhibitory bath for the inhibitory neurons - for the importance of this bath see the SI Appendix of Börgers et al. [1].
+# These currents are specified as distributions using the syntax from [Distributions.jl.](https://juliastats.org/Distributions.jl/stable/starting/) The advantage to this is that a distribution can be given to a call of ``rand()`` and the random number will be drawn from the specified distribution. We'll use this call during the neuron creation step below.
 
 I_base = Normal(0, 0.1) ## base external current for all neurons
 I_driveE = Normal(μ_E, σ_E) ## External current for driven excitatory neurons
@@ -64,15 +71,15 @@ I_bath = -0.7; ## External inhibitory bath for inhibitory neurons - value from p
 # # Creating a network in Neuroblox
 # Creating and running a network of neurons in Neuroblox consists of three steps: defining the neurons, defining the graph of connections between the neurons, and simulating the system represented by the graph.
 
-# ## Define the neurons
+# ### Define the neurons
 # The neurons from Börgers et al. [1] are implemented in Neuroblox as `PINGNeuronExci` and `PINGNeuronInhib`. We can specify their initial current drives and create the neurons as follows:
 
 exci_driven = [PINGNeuronExci(name=Symbol("ED$i"), I_ext=rand(I_driveE) + rand(I_base)) for i in 1:NE_driven] ## In-line loop to create the driven excitatory neurons, named ED1, ED2, etc.
 exci_other  = [PINGNeuronExci(name=Symbol("EO$i"), I_ext=rand(I_base) + rand(I_undriven)) for i in 1:NE_other] ## In-line loop to create the undriven excitatory neurons, named EO1, EO2, etc.
 exci        = [exci_driven; exci_other] ## Concatenate the driven and undriven excitatory neurons into a single vector for convenience
 inhib       = [PINGNeuronInhib(name=Symbol("ID$i"), I_ext=rand(I_driveI) + rand(I_base) + I_bath) for i in 1:NI_driven]; ## In-line loop to create the inhibitory neurons, named ID1, ID2, etc.
 
-# ## Define the graph of network connections
+# ### Define the graph of network connections
 # This portion illustrates how we go about creating a network of neuronal connections.
 
 g = MetaDiGraph() ## Initialize the graph
@@ -92,8 +99,30 @@ for ni1 ∈ inhib
     end
 end
 
-# ## Simulate the network
+#md # !!! note "Alternative graph creation"
+#md #       If you are creating a very large network of neurons, it may be more efficient to add all of the nodes first and then all of the edges via an adjacency matrix.
+#md #       To illustrate this, here is an alternative to the previous block that will initialize the same graph. Just uncomment this code and use it to replace the code from the section above.
+#md #       ```julia
+#md #       g = MetaDiGraph() ## Initialize the graph
+#md #       add_blox!.(Ref(g), [exci; inhib]) ## Add all the neurons to the graph
+#md #       adj = zeros(N_total, N_total) ## Initialize the adjacency matrix
+#md #       for i ∈ 1:NE_driven + NE_other
+#md #           for j ∈ 1:NI_driven
+#md #               adj[i, NE_driven + NE_other + j] = g_EI/N
+#md #               adj[NE_driven + NE_other + j, i] = g_IE/N
+#md #           end
+#md #       end
+#md #       for i ∈ 1:NI_driven
+#md #           for j ∈ 1:NI_driven
+#md #               adj[NE_driven + NE_other + i, NE_driven + NE_other + j] = g_II/N
+#md #           end
+#md #       end
+#md #       create_adjacency_edges!(g, adj)
+#md #       ```
+
+# ### Simulate the network
 # Now that we have the neurons and the graph, we can simulate the network. We use the `system_from_graph` function to create a system of ODEs from the graph and then solve it using the DifferentialEquations.jl package, but for performance scaling reasons we will use the experimental option `graphdynamics=true` which uses a separate compilation backend called [GraphDynamics.jl](https://github.com/Neuroblox/GraphDynamics.jl). The GraphDynamics.jl backend is still experimental, and may not yet support all of the standard Neuroblox features, such as those seen in the Spectral DCM tutorial.
+# We choose to solve this system using the ``Tsit5()`` solver. If you're coming from Matlab, this is a more efficient solver analogous to ``ode45``. It's a good first try for systems that aren't really stiff. If you want to try other solvers, we'd recommend trying with ``Vern7()`` (higher precision but still efficient). If you're **really** interested in solver choices, one of the great things about Julia is the [wide variety of solvers available.](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
 
 tspan = (0.0, 300.0) ## Time span for the simulation - run for 300ms to match the Börgers et al. [1] Figure 1.
 @named sys = system_from_graph(g, graphdynamics=true) ## Use GraphDynamics.jl otherwise this can be a very slow simulation
@@ -115,7 +144,12 @@ fig
 # # Conclusion
 # And there you have it! A complete PING demonstration that reproduces the dynamics of a published paper in a matter of 30 seconds, give or take. Have fun making your own!
 
+# # For further exploration
+# If you want to explore the PING network further, you can try the following:
+# 1. You might notice that the excitatory and inhibitory populations become slightly desynchronized by the end of the simulation, unlike in the original paper. This is because of slight differences in how we implement the excitatory drive and inhibitory bath, which adjusts the overall E/I balance. Try increasing the inhibitory bath or decreasing the percentage of excitatory neurons that receive input and see how this affects the synchrony!
+# 2. This tutorial is written entirely using built-in blocks for the excitatory and inhibitory neurons. If you want to explore the details, try typing ``?PINGNeuronExci`` or ``?PINGNeuronInhib`` in your Julia REPL to see the full details of the blocks. If you really want to dig into the details, type ``@edit PINGNeuronExci()`` to open the source code and see how the equations are written.
+
 # # References
 # 1. Börgers C, Epstein S, Kopell NJ. Gamma oscillations mediate stimulus competition and attentional selection in a cortical network model. Proc Natl Acad Sci U S A. 2008 Nov 18;105(46):18023-8. DOI: [10.1073/pnas.0809511105](https://www.doi.org/10.1073/pnas.0809511105).
-# 2. Traub, RD, Miles, R. Neuronal Networks of the Hippocampus. Cambridge University Press, Cambridge, UK, 1991.
-# 3. Wang, X-J, Buzsáki, G. Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. J. Neurosci., 16:6402–6413, 1996.
+# 2. Traub, RD, Miles, R. Neuronal Networks of the Hippocampus. Cambridge University Press, Cambridge, UK, 1991. DOI: [10.1017/CBO9780511895401](https://www.doi.org/10.1017/CBO9780511895401)
+# 3. Wang, X-J, Buzsáki, G. Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. J. Neurosci., 16:6402–6413, 1996. DOI: [10.1523/JNEUROSCI.16-20-06402.1996](https://www.doi.org/10.1523/JNEUROSCI.16-20-06402.1996)

src/blox/ping_neuron_examples.jl

@@ -1,7 +1,40 @@
 # First, create an abstract neuron that we'll extend to create the neurons for this tutorial.
 abstract type AbstractPINGNeuron <: AbstractNeuronBlox end
 
-# Create an excitatory neuron that inherits from this neuron blox.
+"""
+    PINGNeuronExci(name, namespace, C, g_Na, V_Na, g_K, V_K, g_L, V_L, I_ext, τ_R, τ_D)
+
+    Create an excitatory neuron from Borgers et al. (2008).
+    The formal definition of this blox is:
+
+```math
+\\frac{dV}{dt} = \\frac{1}{C}(-g_{Na}*m_{\\infty}^3*h*(V - V_{Na}) - g_K*n^4*(V - V_K) - g_L*(V - V_L) + I_{ext} + jcn)
+\\m_{\\infty} = \\frac{a_m(V)}{a_m(V) + b_m(V)}
+\\frac{dn}{dt} = a_n(V)*(1 - n) - b_n(V)*n
+\\frac{dh}{dt} = a_h(V)*(1 - h) - b_h(V)*h
+\\frac{ds}{dt} = \\frac{1}{2}*(1 + \\tanh(V/10))*(\\frac{1 - s}{\\tau_R} - \\frac{s}{\\tau_D})
+```
+where ``jcn`` is any input to the blox. Note that this is a modified Hodgkin-Huxley formalism with an additional synaptic accumulation term.
+Synapses are added into the ``jcn`` term by connecting the postsynaptic neuron's voltage to the presynaptic neuron's output:
+```math
+jcn = w*s*(V_E - V)
+```
+where ``w`` is the weight of the synapse and ``V_E`` is the reversal potential of the excitatory synapse.
+
+Inputs:
+- name: Name given to ODESystem object within the blox.
+- namespace: Additional namespace above name if needed for inheritance.
+- C: Membrane capacitance (defaults to 1.0).
+- g_Na: Sodium conductance (defaults to 100.0).
+- V_Na: Sodium reversal potential (defaults to 50.0).
+- g_K: Potassium conductance (defaults to 80.0).
+- V_K: Potassium reversal potential (defaults to -100.0).
+- g_L: Leak conductance (defaults to 0.1).
+- V_L: Leak reversal potential (defaults to -67.0).
+- I_ext: External current (defaults to 0.0).
+- τ_R: Rise time of synaptic conductance (defaults to 0.2).
+- τ_D: Decay time of synaptic conductance (defaults to 2.0).
+"""
 struct PINGNeuronExci <: AbstractPINGNeuron
     params
     output
@@ -43,7 +76,40 @@ struct PINGNeuronExci <: AbstractPINGNeuron
     end
 end
 
-# Create an inhibitory neuron that inherits from the same neuron blox.
+"""
+    PINGNeuronInhib(name, namespace, C, g_Na, V_Na, g_K, V_K, g_L, V_L, I_ext, τ_R, τ_D)
+
+    Create an inhibitory neuron from Borgers et al. (2008).
+    The formal definition of this blox is:
+
+```math
+\\frac{dV}{dt} = \\frac{1}{C}(-g_{Na}*m_{\\infty}^3*h*(V - V_{Na}) - g_K*n^4*(V - V_K) - g_L*(V - V_L) + I_{ext} + jcn)
+\\m_{\\infty} = \\frac{a_m(V)}{a_m(V) + b_m(V)}
+\\frac{dn}{dt} = a_n(V)*(1 - n) - b_n(V)*n
+\\frac{dh}{dt} = a_h(V)*(1 - h) - b_h(V)*h
+\\frac{ds}{dt} = \\frac{1}{2}*(1 + \\tanh(V/10))*(\\frac{1 - s}{\\tau_R} - \\frac{s}{\\tau_D})
+```
+where ``jcn`` is any input to the blox. Note that this is a modified Hodgkin-Huxley formalism with an additional synaptic accumulation term.
+Synapses are added into the ``jcn`` term by connecting the postsynaptic neuron's voltage to the presynaptic neuron's output:
+```math
+jcn = w*s*(V_I - V)
+```
+where ``w`` is the weight of the synapse and ``V_I`` is the reversal potential of the inhibitory synapse.
+
+Inputs:
+- name: Name given to ODESystem object within the blox.
+- namespace: Additional namespace above name if needed for inheritance.
+- C: Membrane capacitance (defaults to 1.0).
+- g_Na: Sodium conductance (defaults to 35.0).
+- V_Na: Sodium reversal potential (defaults to 55.0).
+- g_K: Potassium conductance (defaults to 9.0).
+- V_K: Potassium reversal potential (defaults to -90.0).
+- g_L: Leak conductance (defaults to 0.1).
+- V_L: Leak reversal potential (defaults to -65.0).
+- I_ext: External current (defaults to 0.0).
+- τ_R: Rise time of synaptic conductance (defaults to 0.5).
+- τ_D: Decay time of synaptic conductance (defaults to 10.0).
+"""
 struct PINGNeuronInhib <: AbstractPINGNeuron
     params
     output