Fractal Mode-Locking and Power Law Statistics in Neurodynamics

Chris King – Mathematics
Department, University of Auckland

*Abstract: This short report highlights properties of fractality and self-organized criticality in the Ising model of ferromagnetism and how these ideas can be applied using wavelet transforms to comparisons with the study of self-organized criticality in neurodynamics. Matlab programs are provided to freely replicate the results.*

* *

**Blog Update on the Quantum Picture: Exceptionally Simple? E8, Golden Ratio and Ising Spins**

*
*

**Introduction:** The Ising model [1]
investigates the phase transition between ferromagnetism and paramagnetism
through the Metropolis-Hastings algorithm [2]
run inside a Monte Carlo loop [3].
For negative interaction strengths between spin pairs, an anti-ferromagnetic
spin-glass results in which adjacent spins have lowest energy when their spins
alternate up and down.

Fig 1: The Ising model of
ferromagnetism as an example of phase transition criticality. In the lower figure a cellular
simulation based on a given cell element interacting with the four nearest
neighbours (above and below and to either side) in a rectangular array. The
array is iterated according to the Hamiltonian, in this case for 3000 steps. In
blue magnetization is plotted as a function of interaction strength of
individual spin elements, with the theoretical curve outlined in red. For low
interaction strengths (or equivalently for high thermal energies by comparison
with the interaction strength) we have paramagnetism in which individual
domains of cells of aligned spin are randomly distributed and net magnetization
will occur only in the presence of an external magnetic field. As the
interaction strength increases, there is a critical phase transition to the
ferromagnetic state after which the spins become polarized into large domains
with aligned spins and in the asymptotic limit dominance of either spin up or
spin down results in bifurcation towards a fully saturated state of net
magnetization 1 or -1. At the critical transition value of J = ~0.44
magnetization domains form a fractal distribution similar to states of
self-organized criticality, such as earthquakes and sand-piles, in which the
critical state is maintained by fractal avalanches. In the top row interaction
strength is plotted for both positive and negative interaction strengths on a
larger array that shown in the lower figure. When the interaction strength is negative spins have a lower
energy when adjacent spins are opposite, so we then have the anti-ferromagnetism
of a spin-glass. Negative values also have a critical transition, but this does
not result in net magnetization as the lowest energy states form a
chequer-board having zero net polarization. However, as there are two
complementary chequer-board arrangements, large domains still develop, with
frustration along their boundaries, where the energy cannot be minimized. In
the insets are shown corresponding charts of energy versus magnetization and
interaction strength.

**Methods:** A set of Matlab m-files
performing computational Monte Carlo simulations, developed from [4]
following [5] are provided in the link in appendix 2. The equations supporting the algorithm and some of the theory is displayed in appendix 1.

Generally the computational simulation is performed using only nearest neighbours directly above and below and to either side of an element of a rectangular array, (although more elaborate neighbourhoods are also used in figure 2). This closest neighbour computational process coincides with the Bethe-Peierls approximation
[6]
to Ising spin states in statistical mechanics.

An
array is initiated in a random configuration and is iterated cell by cell,
flipping a set proportion of the spins in each iteration if a random variable
exceeds the energy difference between the cell and its flipped state in
relation to its four neighbours. This provides a thermodynamic model in which
spins will flip to a lower energy state but may also, with exponentially
diminishing probability, become flipped to a higher energy one.

Figures 1 and 2 show the results of this investigation both for the fractal dynamics of the critical phase transition to ferromagnetism with positive interaction energies and for the capacity to form fractal mode-locked states in an anti-ferromagnetic spin-glass state with negative interaction energies. Figure 3 also shows this fractal critical behavior by analyzing the connected regions within the final state of the array.

Fig 2: The effect of increasing
polarization on an anti-ferromagnetic spin-glass in which the lowest energy
state has neighbouring cells having oppositely aligned spins. In the lower
figure (a) is shown the magnetization for a 2-D spin glass in which only the
four closest neighbours interact, as the external field B is increased. There is a strong plateau in which a
chequer-board lowest energy state with adjacent spins opposite is maintained.
(b) There is also a hint of plateuing at a magnetization of ~0.5 in which
¼ of the spins would be aligned one way and ¾ the other. The top
row arrays show states corresponding to various values of the external field.
(c) The theoretical distribution for a 1-D spin-glass with negative interaction
strength (Schr?der) is a ?devil?s staircase? in which mode-locking occurs for
an interval neighbourhood of each rational magnetization state in which a
regular periodic arrangement, such as *uduud*, can result in rational periods
of any order. (d) The 2-D simulation failed to show any further model-locking
even when the interacting neighbourhood was enlarged to include more distant
interactions (with an inverse square 2-D dipole interaction), however in the
corresponding 1-D simulation (with an inverse linear dipole interaction) a
series of plateaus formed at rational magnetization values.

The negative interaction strength regime was also explored using variations on the Matlab algorithm, in both 2-D and 1-D modes, to test for the fractal devil?s staircase of mode-locked states [7] referred to in Schr?der [8]. The lack of additional plateaus was hypothesized to result from energy interaction confinement to nearest neighbours, so variants of the algorithm were developed in which more distant neighbours were also involved in the energy function through dipole interactions in their respective dimensionalities. While the 2-D model failed to demonstrate additional mode-locking states, the 1-D model did then return a series of near stationary intervals surrounding several fractional states of magnetization, consistent with the results in Schr?der.

Fig 3: Log-log plot of
region size versus number of regions of this size range for J = 0, 0.2, 0.4406868, 0.6, 1 coloured blue, green, black,
magenta and red respectively, each for an average of five Monte Carlo runs of
3000 iterations, showing that the critical state has a consistent fractal power
law distribution (black) by contrast with stochastic sub-critical J which lack
large regions and partitioned super-critical J which have a preponderance of a
few large regions, (lower red and magenta sections) containing a small number
of small random incursions (left).

Recent research in neurodynamics by Bullmore?s group [9]
has highlighted the close correspondence between the fractal dynamics of the
critical state of the Ising model, other models such as the Kuramoto model [10]
and self-organized criticality in brain dynamics. We thus extended our use of
the Matlab algorithm to replicate some of the methods used in their work. We
thus adapted a wavelet analysis algorithm to portray the time behaviour of the
Ising model in terms of wavelet amplitudes and the phase correlations between
pairs of channels.

Fig 4: Above: Wavelet transforms
of the summed outputs of an 8x8 sub-array of a 96x96 Ising simulation at the
values J = 0, 0.4406868, 1 for zero external field
demonstrate fractal time variation of a combined channel at the critical value.
The dominance of low frequency wavelet amplitudes at J=1 reflects the
partitioning into large domains. The high frequency contributions at J=0
reflect the small scale randomness. Below: Plots of phase lock of the complex
correlation function between superimposed channels for the same values of J as
above, showing that the fractal nature of the critical state also extends to a
fractal distribution of temporal
phase lock episodes. The black regions for J = 1 result from zero wavelet
amplitudes of one or other channel at higher frequencies causing the argument
of the complex correlation function to be undefined.

In figure 4 is shown an absolute wavelet transform for an
exponential series of frequencies for sub-critical, critical and super-critical
interaction strengths, demonstrating that the fractal nature of the critical
state can be readily detected using wavelet transforms. We defined and used a
complex version of the Morlet wavelet rather than the
Hilbert wavelet transform of Bullmore?s group, but have otherwise followed
their methods, using an Ising simulation on a 96x96 array and then summed the
binary [-1,1] values on each 8x8 sub-array to form 12x12 - effectively 64-bit
channels of ?real? amplitudes – iterated over the last 8192 iterations of
a 12,192 iteration run.

The successive iterates were firstly given a direct wavelet
portrait using the absolute amplitude of the real part of the wavelet
coefficient array to generate the upper series of profiles in figure 4. Pairs of channels were then phase
correlated and the phase angle of the resulting complex argument taken in the
range to generate a
set of two-channel correlation profiles, giving a portrait of the fractal
distribution of the critical state in terms of time intervals of phase
correlation between a pair of channels.

Fig 5: Log-log plot of
number of intervals on which a pair of channels have phase lock with respect to
the length of the phase lock. J = 0, 0.4406868, 1
coloured blue, black, and red respectively. The critical state (black)
approximates a power law distribution, while J = 0 lacks longer phase lock
intervals and J=1 is dominated by large intervals corresponding to the large
domains and very short intervals caused by frustration between competing
domains.

Finally we pooled phase correlation intervals of several pairs of
channels using Bullmore group?s criteria to test for a power law relationship
at criticality over 4 orders of magnitude of interval length (exponentially
double the range of fig 6). These results were not averaged over short time
intervals as in their work, and were performed only for a small number of
pairs, but do show an approximate power law distribution in the log-log plot of
figure 5 as a proof of principle for the Matlab program which can be compared
with their results as shown in figure 6.

Fig 6: (Left): Probability
distribution of phase lock interval (PLI) between pairs of processes at critical
(black line)

and at hot temperature (low
interaction strength) (blue)
plotted on a log-log scale. The black dashed line represents a power law with
slope a~{1:5. (Right) Probability distribution of lability of global
synchronization (the number of channel pairs in phase lock over a given
interval) at critical temperature (black line) and at hot temperature (low
interaction strength) (blue); the black dashed line represents a power law with
slope a~{0:5) (Bullmore et. al.).

**Appendix 1: Summary of the
Theory and Equations**

The Hamiltonian for the
interaction is .

The ratio of probabilities for a
flip is . The spins of particles in which *r* > 1 or *r*
< 1
but greater than a uniformly distributed random number have the potential to be
flipped if a second random number exceeds a given threshold regulating the
proportion flipping in each iteration.

The values for the theoretical
thermodynamic curves shown in figure 1 are outlined below:.

The wavelet coefficients the phase
coherence and the lability
of global synchronization are given by:

**Appendix 2: Matlab Programs
Download Link**

** **

** **

**References**

[4] Gudmundsson J 2008 *Monte
Carlo Method and the Ising Model*

http://www.isv.uu.se/~ingelman/graduate_school/courses/montecarlo/hand-in/jon_emil_gudmundsson.pdf

[5] Fricke T 2006 *Monte
Carlo investigation of the Ising model*

http://web1.pas.rochester.edu/~tobin/notebook/2006/12/27/ising-paper.pdf

[6] Huang K 1963, 1987 __Statistical
Mechanics__, John Wiley and Sons, NY p 321.

[8] Schr?der M. (1993)
Fractals, Chaos and Power Laws ISBN 0-7167-2136-8 p 357.

[9] Kitzbichler M, Smith
M, Christensen S, Bullmore E (2009) *Broadband Criticality of Human Brain
Network Synchronization* PLoS Computational
Biology, 5, e1000314.