In 1952, Hodgkin and Huxley developed a mathematical model to explain the ionic mechanisms underlying the initiation and propagation of action potentials (APs) in the squid giant axon (Hodgkin and Huxley, 1952). In the HH model, each component of an excitable cell is treated as an electrical element. The lipid bilayer is represented as a capacitance (C_m). Voltage-gated ion channels and leak channels are represented by electrical conductance (g_Na, g_K, and g_L denote the maximum conductance per surface area of the sodium, potassium and leak currents). Finally, ion pumps are represented by current sources (I). Explicitly, the dynamics of neurons is described by the following equations:

where a set of four time-dependent variables (V_i, m_i, n_i, h_i) were used to describe the activity of i-th neuron. Here V_i is the membrane potential, m_i and h_i are the activation and inactivation variables of the sodium current, and n_i is the activation variable of the potassium current. V_Na, V_k and V_rest are the corresponding reversal potentials. Typical values of the parameters were chosen as C_m = 1.0 μF/cm², g_Na = 120.0 mS/cm², g_K = 36.0 mS/cm², g_L = 0.3 mS/cm², V_Na = 115.0 mV, V_k = −12.0 mV, and V_rest = 10.6 mV. Iisyn(t) is the total synaptic current, a sum of output currents, Ijout(t-T0), from connected neurons in the network with a synaptic strength w_ij and a signal delay time T₀. Explicitly, the total synaptic current is expressed as Iisyn(t)=∑jwij·Ijout(t-T0). For simplicity we approximated the output current I_j^out(t) as a step function with a duration 0.1 ms and an amplitude Imax·{1+exp[-0.002·Vjpeak(t)]}-1, where V_j^peak(t) is the peak value of the AP of j-th neuron at time t and I_max is the maximum output current from a neuron (Koch, 1999). Typical values of I_max and T₀ used in our simulations are 25 nA/cm² and 9 ms, respectively. The strength or efficacy (w_ij) of synaptic transmission at preexisting synapses is subject to activity-dependent modification, and will be discussed in detail in Section Synaptic Plasticity.

For a small enough time step Δt, the neuronal membrane potential V(t) in Equation (1) can be solved numerically using the Euler method as:

where n(t), m(t), and h(t) are obtained from solving Equations (2) to (4) (Koch, 1999). For a brief derivation, Equations (2)–(4) can be expressed as dzdt=αz(V)(1-z)-βz(V)z = (z∞-z)/τz, where z_∞ = α_z/(α_z+β_z), τ_z = 1/(α_z+β_z), α_z and β_z are voltage-dependent rate constants in Equations (2)–(4), and z represents m, n, or h. When the membrane potential is held at a constant value (such as by voltage clamp), the solution of the three gating equations can be obtained as z(t) = z_∞ + (z₀ − z_∞) exp (−t/τ_z), where z₀ = 0.05, 0.32, and 0.60 respectively for z = m, n, and h. In general, the neuronal membrane potential obtained from Equation (5) is more accurate with a smaller step size Δt. In our experience, Δt = 0.01 ms will be good choices for the size of time step. In this study, for the efficiency of the network simulations and the accuracy of simulation results, we chose Δt = 0.01 ms.

In neural systems, a synapse between two neurons can change its strength in response to either use or disuse of transmission over synaptic pathways (Hughes, 1958). Previously, the Hebbian learning rule was suggested that the synaptic strength could increase if the presynaptic neuron repeatedly and persistently stimulates the postsynaptic neuron to generate APs. More recent experiments have observed a spike-timing-dependent synaptic plasticity (STDP): repeated presynaptic spike arrival a few milliseconds before postsynaptic action potentials leads in many synapse types to long-term potentiation (LTP) of the synapses, whereas repeated spike arrival after postsynaptic spikes leads to long-term depression (LTD) of the same synapse. Previous experiments also demonstrated that postsynaptic APs are initiated in the axon and then propagate back into the dendritic arbor of neocortical pyramidal neurons, evoking an activity dependent dendritic Ca²⁺ influx that could be a signal to induce modifications at the dendritic synapses that were active around the time of AP initiation (Markram et al., 1997). Therefore, the synaptic efficacy can be regulated depending on the precise timing of postsynaptic APs relative to excitatory postsynaptic potentials. The characteristic time intervals for synaptic modifications are found to be 17 ms for facilitation and 34 ms for depression for layer 5 pyramidal neurons in somatosensory cortex (Bi and Poo, 2001). Such a STDP rule was introduced in our simulations by considering a change in the synaptic strength (Δw_ij) due to learning at each time step as Bi and Poo (2001):

where Δτ is the time of the postsynaptic spike minus the time of the presynaptic spike. The parameters τ₊ and τ₋ determine the ranges of pre-to-postsynaptic inter-spike intervals over which synaptic strengthening and weakening occur. A₊ and A₋, which are both positive, determine the maximum amounts of synaptic modifications. Typical values of these parameters in our simulations were A₊ = 0.013, A₋ = 0.005, τ₊ = 10 ms and τ₋ = 9.5 ms.

Alternatively, we also considered the inverse STDP rule that has been observed, for example, in the sensory systems and cerebral cortex of fish (Bell et al., 1999; Zhigulin et al., 2003; Fino et al., 2005). In this case, the change in the synaptic strength due to learning at each time step is expressed as:

Typical values of these parameters in our simulations are A₊ = 0.005, A₋ = 0.013, τ₊ = 9.5 ms and τ₋ = 10 ms.

In this study, the matrix of synaptic strength is asymmetric. A change in the synaptic strength might result from remodeling of synapses in both presynaptic loci and postsynaptic terminals. Saturation of synaptic efficacy can occur after repeated potentiation, and previous experiments have shown that saturation of hippocampal LTP impairs spatial learning (Castro et al., 1989; Moser et al., 1998). Here, for simplicity, we did not introduce saturation for the synaptic efficacy since we studied the simplified case that the learning due to the external signal is turned off before the synaptic efficacy saturates. Nevertheless, saturation could be achieved by an appropriate choice of A₊ and A₋ as a function of the synaptic strength (w_ij), which vanishes at the maximal w_ij. We note that several other stabilization procedures can also terminate learning either when activity levels reach a certain threshold level or invoke a bound on synaptic weight strengths (Nass and Cooper, 1975; Linsker, 1986).

The dynamics of neural networks was studied by considering the HH model (Equations 1–4) in Section The Neuron Model and the synaptic plasticity rules (Equations 6, 7) in Section Synaptic Plasticity. In addition, in a noisy environment, electric noise could play an important role in the neuron dynamics. For example, in the experiment of hippocampal CA3 networks ex vivo by multi-neuron imaging technique (Takahashi et al., 2010), it was observed that neurons can synchronize in a noisy network. The main source of this noise is typically synaptic, resulting from the probabilistic release of synaptic vesicles and bombardment from the myriad of synapses made by other cells. Although neurons' firing frequency rarely exceeds 100 Hz, the combined synaptic activities of a neural network can produce fluctuations on a much faster time scale. Synaptic noise causes abrupt changes in the associated synaptic conductance each time a spike invades the pre-synaptic bouton. Stein's model describes its effect in the evolution of the membrane potential of a given neuron as trains of Dirac delta functions, which are summed up to become Gaussian white noise in the diffusion limit of synaptic input (Stein, 1967). Thus, in the presence of synaptic noise, the membrane potential of a network neuron can be solved using the Euler method as:

where the Gaussian white noise (with the standard deviation D), I_noise, was generated using the Box–Müller transform. The typical value of D in this study is 25 μA/cm².

For the simulation of coupled neural networks, we first built two independent neural layers (layers 1 and 2), each consisting of M = 50 excitatory HH neurons grafted on a square substrate of area 10⁴ (in the unit of soma area). In this study, for the positioning of neurons on each layer, we consider three different distributions, including a random distribution, a grid distribution, and a lognormal distribution. In the case of a lognormal distribution, the x- and y-coordinates of neurons were determined using the function min(100 × f_logn (x|0.5, 0.2), 100), where flogn(x|μ,σ)=(xσ2π)-1exp[-(ln x-μ)2/2σ2] is a lognormal function with the log mean μ and log standard deviation σ, and min(a, b) selects the smaller value from a and b. All neuron pairs have a distance longer than 1 unit length to avoid overlapping of neurons. Our motivation to use the lognormal distribution is to increase the population of neurons with high synaptic connectivity, which serve as hub neurons in a neural network and enhance network-wide synchronicity (Bonifazi et al., 2009).

For each neural layer, based on previous studies (Jia et al., 2004; Lin et al., 2011), we assume that the probability to form a synaptic connection between intra-layer neurons is inversely proportional to the distance between neurons, i.e., pi,jintra=k/ri,j, where r_ij is the distance between neurons i and j. Such a connection probability is valid if the process for a neuron to find another one is by a random search on a two-dimensional plane. The coefficient k in general depends on the level of activity in the network, such as local concentration of neurotrophins (Vicario-Abejón et al., 1998), and we used k = 0.005 in our simulations. The total number of intra-layer connections on a neural layer is denoted as N_c(N_c1 for layer 1, and N_c2 for layer 2).

For the coupling between two layers, as shown in Figure 1, we considered two types of inter-layer connections, including a random connection and a preferential connection. In the former case, the inter-layer synapses were established randomly between neurons with a fixed probability; while in the latter case, the probability to form an inter-layer synapse was proportional to the number of existing intra-layer connections of neurons i and i′. In other words, we assumed pi,i′inter=cici′/(M−1)2, where M is the number of neurons per layer and c_i is the number of intra-layer connections of neuron i (neurons i and i′ belong to different layers). We differentiated between the connection from neuron i to neuron j and that from neuron j to neuron i, i.e., all intra- and inter-layer connections are directed. In our simulation, we assumed an equal number (N_i) of inter-layer connections in either direction.

To simulate the dynamics of the coupled neural network, initially we set the membrane potential of neurons to have a Gaussian distribution with a zero average and a standard deviation of 5 mV and the synaptic strength between neurons to have a Gaussian distribution with an average of 0.025 and a standard deviation of 0.01. The simulation of network activities was divided into two parts, each containing a learning phase for 2 s and a recall phase for 3 s. The synaptic strength of each connection was updated according to learning rules (Equations 6 or 7) in the learning phase, but remained a constant value in the recall phase. In the first part of simulation, we studied the dynamics of two independent layers by varying the number of intra-layer synaptic connections (N_i) and the maximum amounts of synaptic modification (A₊ and A₋), and obtained the phase diagram of each independent layer. In the second part of simulation, we investigated the interaction of these two layers by varying the number of inter-layer synaptic connections. In addition, the inter-layer connections were disconnected at the end of the recall phase in the second part of simulation, and we simulated the firing dynamics of the network for another 3 s (without learning) to investigate the stability of the induced synchronous firing state (SFS). We note that all parameters in our simulations were assumed to be their typical value unless otherwise specified. For each set of parameters, we simulated the network activities 30 times by using different seeds to measure the average value of related physical properties.

The hypothesis that the brain computes information using neural synchronization has been supported by a mounting number of experimental evidence (Rodriguez et al., 1999; Fries, 2009). In this study, we first studied the noise-driven synchronous firing dynamics of a developing neural network (a single layer of 50 HH neurons) and derived its phase diagram. To begin with, there was no connection among neurons in the network. As the network developed with culturing time, intra-layer connections were established between neurons with the probability pi,jintra. As demonstrated in our previous study (Lin et al., 2011), as network connectivity increased, more and more neurons became active due to the presence of noise. Two neurons (i and j) were considered as synchronized in their firing pattern if both of them are active (fired at least once in the learning phase) and their time correlation (TC_ij ≡ cov(V_i(t), V_j(t))/σ_iσ_j, where cov means covariance and σ_i is the standard deviation of V_i) is greater than 0.2 (Chao and Chen, 2005; Lin et al., 2011). The selected threshold value of time correlation for neural synchronization is usually not high since synchronous neuron firing only requires that the membrane potentials of synchronous neurons reach the threshold potential at the same time. By inspecting our simulation results, we selected 0.2 (which is not rigid) as the threshold value for indicating synchronous firing of neurons. Although there are some false positives in identifying synchronous neurons by using a low threshold value, they can be avoided by adding a constraint in the root mean square deviation of neurons' firing time to ensure that neurons' firing events are within a limited time span. Moreover, we defined an order parameter Ψ_s for the synchronization of neural activities as the number of synchronized neuron pairs divided by the total number of connections. In the present study, a network state is considered as a SFS if the time average of Ψ_s (〈Ψ_s〉) in the recall phase is greater than 0.95, a transition state (TS) if 0.4 ≤ 〈Ψ_s〉 ≤ 0.95, and a background activity state (BAS) if 〈Ψ_s〉 < 0.4.

In non-linear dynamical systems, it is known that non-trivial effects of noise could lead to synchronization of the system (Neiman et al., 1998; Zhou and Kurths, 2003). Such noise-driven synchronization is a topic of relevance to neuroscience, and has been experimentally observed in animal neocortical neurons (Mainen and Sejnowski, 1995). By varying the degree of network development (i.e., the degree of network connections, N_c), here we investigate the effect of a Gaussian noise (of standard deviation D) on the synchronous firing of a neural network with STDP learning rule, which consists of 50 neurons randomly distributed on a square substrate. In our simulations, no synchronous firing was observed for networks at early developmental stages (small N_c). As the network further developed, synchronous firing was observed for N_c greater than a threshold value (Nc*). The value of Nc* was found to increase as D decreases, and no synchronous firing was observed for D = 18.5. Figure 2 delineates the dependence of Nc* on the standard deviation (D) of I_noise, in which the dotted curve is a fit of data points using Nc* = 2000/(D − 18.5)^0.4. Our results suggest that noise of a large D value can effectively enhance the coherence of the spike trains in a network of coupled neurons.

For the above neural network, its phase diagram is shown in Figure 3 for A₋ = 0.005 and D = 25 μA/cm². The data points in Figure 3 were obtained from calculating the synchronization order parameter (as shown in the inset) and the dotted lines were fitted curves for these data points. In general, for a neural network with the STDP rule, SFS is usually found in the region of large A₊ and large N_c, while BAS exists in the small A₊ and small N_c region of the phase diagram. This observation is consistent with the results in our previous study (Lin et al., 2011). Recent experiments on the synaptic plasticity in mouse barrel cortex have shown that the learning efficacy declines and disappears with age (Banerjee et al., 2009; Kaczorowski and Disterhoft, 2009). Therefore, the BAS network might be associated with neural systems that are aged (small A₊) or undeveloped (small N_c). The inset of Figure 3 shows the value of 〈Ψ_s〉 as a function of N_c (intra-layer connection), which increases sharply for N_c > 700 and saturates around N_c = 1000. Such a narrow TS region (about a range of 300 in N_c) sandwiched by BAS and SFS is also seen for other sets of parameters, as indicated by Figure 3. To further discuss the finite size effect on the transition from BAS to SFS, in Figure 4, we calculated the average synchronization order parameter as a function of the average number of connections per neuron (N_c/M) for different network sizes (M = 50, 80, 100, and 120). In general, data points in Figure 4 collapse into a single curve, suggesting a small finite size effect on the BAS-SFS transition. However, a slightly sharper transition is observed for networks of larger sizes. In our model, SFS is observed if the average number of connections per neuron (N_c/M) is greater than 20 for networks of size between 50 and 120 neurons.

To further our investigation, we considered the firing dynamics of a network consisting of two weakly coupled neural layers. Each layer is composed of 50 neurons randomly positioned on the substrate. There are N_c1 intra-layer connections for layer 1 and N_c2 intra-layer connections for layer 2. Depending on the number of intra-layer connections, these two layers could be in a SFS, a TS, or a BAS before we coupled them together, as described in Section Noise-Driven Synchronous Firing of a Developing Neural Network. Here we first studied the interaction of a SFS layer (layer 1, N_c1 = 1000 and A₊ = 0.013) and a BAS layer (layer 2, N_c2 = 1000 and A₊ = 0.010), which were randomly connected with N_i inter-layer connections. The inter-layer connections were introduced at t = 0 and the learning phase lasted for 2 s, followed by a recall phase for 3 s. For each N_i, 30 realizations were simulated with different seeds to calculate the synchronization order parameter of layer 2. For each realization, the inter-layer connections were varied, while the intra-layer connections were fixed. To study the robustness of the induced synchronization of the BAS layer by a small inter-layer coupling, we then removed all inter-layer connections from this network, and simulated the activities of the network for 3 s (in a recall phase). The average synchronization order parameter of layer 2 is displayed in Figure 5 for various numbers of inter-layer connections, N_i. For a coupled network with N_i ≤ 120, layer 2 did not reach a SFS after 2 s learning, and the order parameter dropped significantly (less than 0.6) after the removal of all inter-layer connections. For N_i = 140, layer 2 has reached a SFS after 2 s learning, but the order parameter fluctuated largely after removing inter-layer connections. For N_i > 160, after 2 s learning, layer 2 reached a SFS and fired synchronously even after the inter-layer coupling was disconnected. It is known that both SFS and BAS are important states for neural computations. Our simulations confirm that SFS occurs for networks with large connectivity or large learning efficiency. These SFS networks can also be desynchronized by various control mechanisms, such as the inhibition control using GABAergic neurons (Treviño, 2016). On the other hand, for those network structures that generate BAS, there is no known efficient mechanism to induce synchronization in these networks without changing the network connectivity or inter-neuron interactions. Therefore, we consider region for SFS in the phase diagram as a better substrate of neural computation than that for BAS. This investigation demonstrates a possible mechanism for the repair or an enhanced learning of a BAS neural layer (possibly resulting from aging, immaturity, or external damage) by coupling it with a SFS layer. Furthermore, we investigated the induced synchronization of a BAS layer by its coupling to a SFS layer for the case of N_i ≥ 180 in Figure 5. The BAS layer was considered as an aged neural layer, which had a learning efficacy (A₊ = 0.010) smaller than that (A₊ = 0.013) of a normal SFS layer with the same degree of intra-layer connectivity (N_c1 = N_c2 = 1000). In Figure 6, we showed the time series of neuron firing and the average neuronal membrane potential (<V₂>) of layer 2 before coupling (a), after coupling (b), and after disconnecting the coupling (c), all of which were recorded in the recall phase. Before coupling, layer 2 showed a random firing pattern and its average membrane potential was usually smaller than the firing threshold (≅ 6.9 mV). After coupling, layer 2 was driven to fire synchronously and its average membrane potential showed a periodic spiking pattern. After the removal of all inter-layer synaptic connections, this periodic spiking persisted for a long time (as shown in Figure 5). We note that, as described earlier for Figure 5, there is no learning phase after decoupling the inter-layer connections.

In Figure 7, we showed the phase diagram of layer 2 in a coupled network of two layers, in which N_c1 = 1000 and N_c2 is variable. Before coupling, layer 1 was in the SFS, layer 2 was in the BAS, and both layers had a learning efficacy A₊ = 0.013. By varying the number of inter-layer connections, N_i, layer 2 was driven to become a SFS layer after learning. The dotted lines showed phase boundaries in the phase diagram that were linearly fitted using the data points in Figure 7. For N_c2 = 100, the value of N_i is about 50 for a transition from a BAS layer to a TS layer, and is about 220 for a transition from a TS layer to a SFS layer. These critical values of N_i decrease linearly with N_c2. It is also observed that the range of N_i for the existence of a TS layer becomes narrower at larger N_c2.

To compare the learning effects of the STDP and the inverse STDP in inducing network synchronization, as shown in Figure 8, we calculated the synchronization order parameter of layer 2 using these two learning rules in a single layer network (a) and a network consisting of two coupled neural layers (b). For the case of a single layer network with various values of N_c, as shown in Figure 8A, we found no difference in the effects of the STDP or the inverse STDP rules since there were roughly equal amount of firing events with positive Δτ or negative Δτ (defined in Equations 6, 7), and thus their contribution to the change of average network synaptic weight was similar with these two learning rules. For the case of two coupled layers in which layer 1 was a SFS layer and layer 2 was a BAS layer, there were more synaptic currents from the direction of layer 1 to layer 2 than the opposite direction. Such an asymmetry led to larger enhancements in the average inter-layer synaptic weight with the STDP learning rule, as shown in the inset of Figure 8B. For N_i > 200, the synchronization order parameter of layer 2 saturated with either learning rule. Nevertheless, there still was a large difference in the average inter-layer synaptic weight of the network with these two learning rules.

In Figure 9, we investigated the interaction of two neural layers (N_c1 = N_c2 = 1000, and A₊ = 0.013 for both layers) having the same frequency but different phases in their synchronous firing pattern. Before the coupling, both layers were in the SFS, and the synchronous firing pattern of layer 1 was lagging behind that of layer 2 by 0.005 s (or a phase difference of 2.9 radians). The inter-layer connections were introduced at t = 0. The learning phase lasted for 2 s, which was followed by a recall phase for 3 s. After introducing the inter-layer coupling, in the recall phase, we observed a phase shift of layer 1 due to the synaptic currents from layer 2 as shown in Figure 9. We note that the data points in Figure 9 were obtained from 5 realizations with the same initial phase difference between the two layers, and the dotted line is a linear fit of data points. As N_i increases, the phase difference between the two layers diminishes. At N_i > 50, the two layers almost fired simultaneously.

As neural synchronization plays an important role in brain's information processing, we would like to explore some possibilities in the positioning of neurons on each layer and in the mechanism of inter-layer connections by designing coupled neural networks, and find out designs that may strengthen their synchronization activities. The realization of these designs might rely on experimental controls of neurotrophins (Vicario-Abejón et al., 1998) as well as neuron positioning related genes, such as the reeler gene and the mouse disabled1 gene (Stanfield and Cowan, 1979; Howell et al., 1997). Similar to those networks discussed in Section Interaction of Two Neural Layers, the coupled networks under our investigation here consist of two coupled neural layers, each of which is composed of 50 neurons. For each layer, as shown in Figure 10, we considered three types of neuron positioning, including a random distribution (a), a grid distribution (b), and a log-normal distribution (c). For these three types of neuron positioning, the degree of intra-layer synaptic connections of each neuron was calculated using the connection rule described in Section Simulating Coupled Neural Networks and averaged over 1000 seeds, as shown in Figure 10D. The total number of intra-layer connections is 1000 for all three cases. It is seen that the distribution of intra-layer connections is almost the same for the cases of (a) and (b), and we expect little difference in their synchronization behavior. On the other hand, for the lognormal positioning of neurons, there is a larger population of neurons which have a large intra-layer connectivity, and we expected to see some enhancement in the synchronization behavior of the network.

Next, we studied the network synchronization behavior by considering two different mechanisms of inter-layer connections as described in Section Simulating Coupled Neural Networks: (a) the random connection and (b) the preferential connection. For the coupled network (N_c1 = 1000, N_c2 = 300) with a random positioning of neurons and N_i inter-layer connections, in Figure 11, we compared the synchronization order parameter of layer 2 for these two connection mechanisms. After the two layers were momentarily coupled, it was found that the value of <Ψ_s2> increased due to the interaction between the two layers. For 100 < N_i < 200, the enhancement in the synchronization order parameter of layer 2 was larger in the case of the preferential connection. Furthermore, we considered six designs of coupled neural networks (three positioning schemes × two connection schemes) and investigated the synchronization enhancement of a BAS layer by its momentary coupling to a SFS layer. As demonstrated in Figure 12, the largest enhancement was observed for the scenario with a lognormal positioning of neurons and the preferential inter-layer connections, and the smallest enhancement was observed for the scenario with a grid positioning of neurons and the random inter-layer connections. These results are valid for networks with the STDP or the inverse STDP learning rules.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.