Ordering Dynamics in Neuron Activity Pattern Model: An insight to Brain Functionality

2.1 MB·English

by  Awaneesh Singh

Ordering Dynamics in Neuron Activity Pattern Model: An insight to Brain
Functionality
Awaneesh Singh1, Jasleen Gundh2, and R.K. Brojen Singh2
∗
1School of Physical Sciences, Jawaharlal Nehru University, New Delhi−110067, India.
2School of Computational and Integrative Sciences,
Jawaharlal Nehru University, New Delhi-110067, India.
We study the ordering kinetics in d = 2 ferromagnets which corresponds to populated neuron
activities with long-ranged interactions, V(r) ∼ r−n associated with short-ranged interaction. We
present the results from comprehensive Monte Carlo (MC) simulations for the nonconserved Ising
modelwithn≥2. Ourresultsoflong-rangedneuronkineticsareconsistentwiththesamedynamical
5
behavior of short-ranged case (n > 4). The calculated characteristic length scale in long-ranged
1 interaction is found to be n dependent (L(t) ∼ t1/(n−2)), whereas short-ranged interaction follows
0
2 L(t) ∼ t1/2 law and approximately preserve universality in domain kinetics. Further, we did the
comparativestudyofphaseorderingnearthecriticaltemperaturewhichfollowsdifferentbehaviours
n
of domain ordering near and far critical temperature but follows universal scaling law.
a
J

2 INTRODUCTION optimize information storage in the system [14] and dy-
2
namic range of the system to response the input signal

[15]. Thecoreofthepaperfocusesonlong-rangeinterac-
] The spiking activity in complex neuron network in
C tions which dynamically explains the rapid movementof
brain is dynamic (far from equilibrium) [1], exhibit
N the signal information inside the brain. The fast emer-
nonequilibrium critical dynamics [2] and the criticality
. gence of long-range interactions, mimic the rapid neu-
o in it has been used to characterizebrain signals [3]. The
ronal interactions in the brain. Further, we study the
i reason could be when such system is quenched from a
b phaseorderingdynamicsinneuronswithaspecificinter-
- homogeneous phase to a broken-symmetry phase, it be-
est to understand the role of the range of inter-neuron
q comesthermodynamicallyunstable. Thesubsequentfar-
[ interactions. We addresstwoimportantquestionsinthis
from-equilibrium evolutionofthe systemis characterized
context via kinetic MC simulations: (a) What is the
1 by the emergence and growth of domains enriched in
growthlawfororderingphasesofneurons? Isthegrowth
v the new equilibrium phases. This nonequilibrium evolu-
5 lawindependentofthe rangeofinteraction? (b)Whatis
tion, usually called kinetics of phase ordering or domain
7 the morphology for ordering phases of neurons, as mea-
growth,hasbeenthesubjectofmuchactiveinvestigation
5 sured by the correlation function and structure factor?
5 [4]. The domain morphology is quantified by the time
Is it comparable for all interaction range? We will be
0 dependence ofthe domainscale L(t), wheret is the time
providing the answers to the above questions from our
. after the quench. There is a good understanding of do-
1 extensive MC simulations.
0 main growth kinetics in pure and isotropic systems with
5 short-rangedinteractions,where the domain scale shows
1 a power-lawbehavior, L(t)∼tφ [5, 6]. For the case with
:
v nonconserved order parameter, e.g., ordering of a ferro-
i magnetintoupanddownphases(spin-flipGlauber-Ising
X
NEURON ACTIVITY PATTERN MODEL
model [7]), one has φ = 1/2 [8, 9]. On the other hand,
r
a for the case with conserved order parameter, e.g., phase
separation of a binary (AB) mixture into A- and B-rich
Brain can be considered as a complex network of
domains (spin-exchange Kawasaki-Ising model [10]), we
neurons which can be mapped onto ferromagnetic Ising
haveφ=1/3whengrowthisdrivenbydiffusion[11,12].
model with two spin interactions [16], where random fir-
Apart from the domain growth laws, experimentalists
ing or non-firing of neuron can be represented by two
arealsointerestedinquantitativefeaturesofthe domain
states of a spin, s=+1 for firing and s=−1 for rest or
morphologies. Animportantexperimentalquantityisthe
non-firingneurons[17]. Eventhoughneuronactivitypat-
time-dependentcorrelationfunctionC(~r,t)oritsFourier
tern model is far from equilibrium dynamic model [18],
transform,structurefactorS(~k,t)(~k beingthewavevec-
Ising model canbe serveanexcellentmodel to deal with
tor [4, 5]. Most of the studies are concentrated around
critical phenomena of neuron activity pattern [1]. The
the nearest-neighbor(short-range)inter-molecular inter-
largenumberoflocal(shortrange)interactionofneurons
action.
[19] and significant amount of global (long range) inter-
The neuron activities in brain at critical point are be- action of neurons [20] are main basis of neurons commu-
lievedtobeeffectiveforthelongdistancecommunication nicationinbrainnetwork[21]. Weconsiderthefollowing
ofthe neurons[13]becauseofcouplingandvariabilityto Hamiltonian of two dimensional Ising system which in-
2
corporates long-rangedspin (neuron) model (LSM), √N
≈ lim J2d(=2) drg(r)r3−n, (3)
H =− J(rij,n)sisj, si =±1, (1) N→∞ Z1
<Xij> where, g(r) is the pair distribution function such that
where J is the coupling strength, n characterizes the g(r)≈1 for r≫1. Then the equation (3) becomes,
range of the interaction, r = |~r −~r |, and s denotes
ij i j i
2ln(N) for n=4,
the spin variable at site i. We consider two state spins:
dsien=ot+es1adednoowtnes-sapninu(pp-aspssinive(ancetuivreonnse)u.roWnes)caonndsidseir=on−ly1 U(n)≈Nl→im∞J n1−41n4/(cid:0)41N−2−Nn/22−n/2(cid:1) ffoorr 0n≤>n4,<4(4)
−
a ferromagnetic case, where J > 0 always. The case 
From equation (4), one can see that U(n) is finite for
where J can be both > 0 (ferromagnetic) and < 0 (an-
n ≥ 4 when N → ∞, and the asymptotic behaviour of
tiferromagnetic)is relevantto spin glasses. We associate
finite critical temperature T [23],
stochastic dynamics with the Ising model by placing it c
in contact with a heat bath. The appropriate dynam- J J 4
T (n)≈ U(n)≈ , (5)
ics for the phase ordering problemis spin-flip kinetics or c
k k (cid:18)n−4(cid:19)
B B
Glauber kinetics.
where, k is Boltzmann constant. This shows that
B
T (n)∝1/n for short-ranged potential (n>4), whereas
c
for long-ranged potential, T depends on the size of the
c
system N as well as n given by,
J 4
T (n,N)≈ N2 n/2, (6)
c −
k (cid:18)4−n(cid:19)
B
andT divergeswith system size. Similarly,one canalso
c
calculate other thermodynamical parameters such as in-
ternal energy, entropy, free energy per particle (neuron)
etc. at this asymptotic limit.
DETAILS OF SIMULATION
Since it is very difficult to obtain exact analytical so-
lution of this problem, we straightforward implement a
MC simulation of the Ising model with spin-flip kinetics
FIG. 1. Plot of < s > vs. T for n= 2, 3, 4, 6, and 12 as to understand the behaviour. In a single step of MC dy-
i
indicated. Themagnetization drops-offsharplynearthecrit- namics, we choose a spin at random in the lattice of dis-
icaltemperature(Tc)andthenvanishesto0inthedisordered tributionof spins. The change inenergy∆H that would
high-temperaturephase.
occurifthespinwasflippediscomputedwiththestepof
acceptance or rejection based on Metropolis acceptance
If we consider interacting neurons through slowly de- probability [27, 28] given by,
caypotentials[22]. Inordertocapturethermodynamical
parameters, one can define the following potential func- exp(−β∆H) if ∆H ≥0,
P = (7)
tion (obeying power law functional form) [23, 24], (cid:26) 1 if ∆H <0.
U(n)= lim 1 N J . (2) where, β = (kBT)−1 denotes the inverse temperature.
N→∞N i,Xj;i=j rinj One Monte Carlo step (MCS) is completed when this
6 algorithm is performed N times (where N is the total
Here n=d+σ =2+σ [25] for two dimensional system. number of spins), regardless of whether the move is ac-
Forshort-rangedinteraction,σ >2andthesystemsizeis cepted or rejected. All our simulations have been per-
notmuchimportant,whereasforlong-rangedinteraction, formedonad=2latticeofsizeL2 (L =512)withperi-
s s
0<σ <2 and it depends on the system size. Since, the odic boundary conditions in both directions. The statis-
size ofthe systemis N,rescalingJ →J/N to the Curie- ticalquantitiespresentedhere(e.g.,correlationfunction,
Weissmodel[24]andusingEuler-McLaurinsumformula structure factor) are obtained as averages over 10 inde-
[26] for N ≫1, equation (2) can be written as, pendent runs. Each run starts with a randomly-mixed
√N √N state with equal numbers of up (fired) and down (inac-
1
U(n)= lim J , tive)spins(neurons),whichcorrespondstoa meanmag-
N→∞ Xx=1Xy=1(x2+y2)(n−2)/2 netization m=hsii=0.
3
vanishes, whereas below T it takes a nonzero value, in-
c
ducing the typicalbehavior ofa ferromagnet. Therefore,
the physical properties of such systems and so its phase
statesdependonthevalueforthemagnetization,thepa-
rameter which is termed as order parameter: an ordered
phaseinwhichthespinsarealignedappearswhenm6=0,
while m = 0 implies a disordered (or symmetric) phase.
Since, T ’s for n ≥ 4 are very close to each-other and
c
hence, exhibit qualitatively similar behavior (explained
shortly). We thus consider n < 4 cases for the long-
ranged interaction. For each value of n, we cut-off the
interaction at r = (2.5)6/n to accelerate our simulation
c
[29]. We stressthatthe simulationsarenumericallyvery
demanding for larger cut-offs. We compute several sta-
tistical quantities to characterize the system. These are
described as follows.
FIG.4. (a)PlotofC(r,t)vs. r/Latt=500forn=2,3,4,6.
(b)Plot of S(k,t)L−2 vs. kL, corresponding to thedata sets
FIG. 2. Evolution snapshots of domain coarsening for n=2, in (a). The reasonably good data collapse shows that the
3,and6quenchedatT =1belowTcattimet=100,500. The scalingfunctionsdonotdependontheinteractionrange. The
snapshotsare obtained from aMonteCarlo (MC) simulation solid line denotes the OJK function in Eq. (12) [33].
of ordering kinetics in ferromagnetic system. The details of
theMC simulation are provided in thetext.
The domain coarsening is characterized by a growing
time-dependent length scale L(t). The domain morphol-
ogydoesnotchangewithtime,apartfromascalefactor.
A direct consequence of the existence of a unique length
scale is that the system exhibits a dynamical-scaling in
the correlation function and structure factor. We com-
pute the time-dependent correlation function:
C(~r ,~r ;t)≡hs s i−hs ihs i. (8)
i j i j i j
Here, the angular brackets denote an averaging over the
independent initial ensemble and different noise realiza-
FIG. 3. (a) Scaling plot of C(r,t) vs. r/L for a phase or- tions. As the systemistranslationallyinvariant,the cor-
dering dynamics in d = 2 for n = 2. The data sets (for relation function depends only on~r =~r −~r :
j i
t = 100,200,500) collapse onto a single master curve. (b)
Similar plot for n=6. C(~r ,~r ;t)=C(~r ,~r +~r;t)=C(~r,t). (9)
i j i i
Usually most experiments study the structure factor,
We study LSM for several values of n, namely 2, 3, which is the Fourier transform of the real-space corre-
4, 6, and 12. The critical ordering temperatures, T (n) lation function:
c
for each n case is shown in Fig. 1, where the charac-
teristic behavior of spontaneous magnetization (< si >) S(cid:16)~k,t(cid:17)=Z d~r C(~r,t)ei~k·~r. (10)
is plotted against temperature (T). As expected, T (n)
c
(dottedlines)increaseswithdecreasingnasevidentfrom Since the system is isotropic, we can improve statistics
equation (5). Above T the spontaneous magnetization bysphericallyaveragingthe correlationfunction andthe
c
4
structure factor. The corresponding quantities are de- with each other and there is no spatial structure in the
noted as C(r,t) and S(k,t), respectively. The correla- evolution morphology. For larger values of n, we see the
tion function and structure factor obey the dynamical emergence and growthof domains of up-spin (marked in
scaling forms: black)anddown-spin(unmarked). These domains inter-
act and annihilate, resulting in coarsening of the char-
C(r,t) =g[r/L(t)],
acteristic length scale, and therefore, domain patterns
S(k,t)=L(t)df[kL(t)]. (11) at different times look statistically similar, apart from a
global change of scale. The domain size at a fixed time
Here, g(x) and f(p) are scaling functions; r is the sep-
(e.g., t=500) is smaller for larger values of n.
aration between two spatial points; k is the magnitude
InFig.3,weshowascalingplotofthecorrelationfunc-
of the wave vector; and d is the system dimensionality.
tion, defined in Eq. (8). We plot C(r,t) as a function of
ThecharacteristicdomainsizeL(t)isobtainedasthedis-
thescaleddistancer/Latthreetimes,asindicated. Fig-
tanceoverwhichthe correlationfunctiondecaysto some
ure 3(a) corresponds to n=2, and Fig. 3(b) shows data
fraction (say half) of its maximum value [C(r,t) = 1
for n = 6. The data sets show a neat scaling collapse,
at r = 0]. There are several other suitable definitions
confirming the validity of dynamical scaling.
for computing L(t), e.g., first zero-crossingof C(r,t), in-
verseofthefirstmomentofS(k,t). Inthescalingregime,
allthesedefinitionsdifferonlybyconstantmultiplicative
factors [30–32].
FIG. 5. Time-dependence of the characteristic length scale
L(t) for n = 2,3,4,6, plotted on a log-log scale. The lines
ofslope1/2indicatesthepower-lawgrowthregimesexpected
for phase ordering in d=2 ferromagnetic system.
FIG. 6. Evolution snapshots of phase ordering systems in
d=2forn=2,3,and6. ThesystemisquenchedatT ≃T .
NUMERICAL RESULTS c
The MC simulations are described in thetext.
In Fig. 2, we show the evolution snapshots obtained Let us next discuss whether the evolutionmorphology
from our MC simulations for n = 2, 3, 6 with T = 1 depends on the range of the interactioncharacterizedby
(< T , see Fig. 1) at t = 100, 500 MCS. At low tem- n. Fig.4 shows a comparisonof the scalingfunctions for
c
peratures,energetic effects are dominantand the system four different n values (n=2, 3, 4, 6) at a time t=500,
minimizesitsenergybyorderingthespinsparalleltoeach when the system is already in the scaling regime. In
other. In the absence of an external field (e.g., magnetic Fig. 4(a), we plot the scaled correlation functions. The
field, h = 0), the activated neuron (up-spin) and inacti- reasonably good data collapse suggests that the scaling
vated neuron (down-spin) states are equivalent. In the functions do not depend on the interaction range. The
mean-field (MF) limit, i.e., n = 0, all the spins interact solidline inFig.4(a)denotestheanalyticalresultdueto
5
Ohta et al. (OJK) [33], who studied ordering dynamics or zero order parameter states, the system near T is a
c
in a ferromagnet. The OJK function is kindoffractalblendofallthree[40]. Howeverthecluster
C(r,t)= 2 sin−1(e−r2/L2). (12) tshizeermisallaregneerrgfoyr(skmaTll)eronf.thReecsaylslttehmatisatloTw,=t1hu(s≪spTicn)s,
π B
try toobtainminimal energyby forming domainswith a
(The corresponding result for the case with vector or-
global ordering: m=+1 or -1.
derparameterhasbeenobtainedbyBrayandPuri[34].)
Our correlation-function data is in excellent agreement
with the OJK function, showing that the phase ordering
dynamics for n < 4 lie in the same dynamical univer-
sality class as that for n > 4. In Fig. 4(b), we plot the
scaled structure factor [L 2S(k,t) vs. kL] for the same
−
time as in Fig. 4(a). Again,the data sets collapse neatly
onto a single master curve, confirming the scaling form
in Eq. (11). The scaling function is in excellent agree-
ment with the corresponding OJK function. Notice that
the structure factor, for large values of k, follows the
well-knownPorod’s law, S(k,t)∼k−(d+1), which results FIG. 7. (a) Plot of C(r,t) vs. r for n = 2, 3, 6 at t = 500
from scattering off sharp interfaces [35, 36]. The scaled MCS.(b)SaclingplotofC(r,t)vs. r/L,correspondingtothe
correlation function and structure factor, in congruence data sets in (a). The solid line denotes theOJK function for
withscalefreebehaviouroffunctionalbrainnetworks[37] thescaled correlation function in Eq. (11).
depicts the universality of the interaction mechanism in
both short and long range interactions in brain. Finally, in Fig. 7(a), we show the plot of correlation
InFig.5,weturnourattentiontothetime-dependence function [C(r,t) vs. r] corresponding to the evolution
of the domain size. We plot L(t) vs. t on a log-log scale shown in Fig. 7 at t = 500. Note that the decay range
for n = 2, 3, 4 and 6. Here, the data sets are consistent of the correlation function is larger for smaller n. Fig-
withthe Cahn-Allengrowthlaw,L(t)∼t1/2−thereisno ure 7(b) shows the scaling plot of the data sets in (a). A
sign of a crossover in the growth law at n = 4, as pre- reasonable data collapse confirms the dynamical scaling
dicted by Bray [38]. Bray has used the renormalization and clarifies that the system for each n belongs to the
group (RG) approach to study ordering dynamics with same dynamical universality class.
long-rangedinteractionsoftheformrd+σ with0<σ <2.
In our case, d = 2 and σ = n−2. Bray argues that the
SUMMARY AND DISCUSSION
long-ranged interactions are relevant for 0 < σ < 2 or
2 < n < 4, and irrelevant for n > 4. The corresponding
growth law is Letusconcludethispaperwithasummaryanddiscus-
sionof the results presentedhere. We study the effect of
t1/(n 2) for 2<n<4,
L(t)∼ − interaction range on the morphology of neuron activity
(cid:26)t1/2 for n>4,
pattern. Thelong-rangedandshort-rangedinteractionof
withpossiblelogarithmiccorrections. Aswecanseethat neurons could be the main basis of how brain performs
our numerical results are not consistent with this pre- complicated functions at fundamental level. Neuron ac-
diction. The only difference as n is varied is that we tivities as well as wiring and rewiring of the neurons in
have faster growth(higher prefactors)for smaller n, cor- the network subjected to heat bath depend on the range
responding to more long-ranged interactions. The fast of interaction which are reflected in the dependence of
dynamics of long-range interactions signifies the path of criticaltemperatureT andmagnetizationonn. Thedo-
c
information processing and neuronal connections in the main sizes of the neurons in short range interaction at
brain. The longer persistance of long-rangedneural con- far-lower critical temperature are smaller; some are iso-
nectionscouldgivesensetoclusteringbehaviourofneural lated and numbers are more as comparedto long ranged
circuitry , specifically during learning of a specific task, interaction for any time domain. However, the domain
new synaptic connections tend to form in vicinity of old dynamics both in short and long ranged interaction sys-
connections related to that task [39] making it more ro- tem is quite different as comparedto far criticaltemper-
bust. aturedynamicsdue toemergenceofmorerandomnessin
InFig.6,weshowtheevolutionoftheorderparameter thedomainorganizationinthesystem. Thisleadstothe
(m)nearcriticaltemperature(T ≃T )fromadisordered changeindomaingrowthlawsoftheneuronsinshortand
c
initial state (m=0) for n=2, 3, and 6 respectively. At longrangedinteractioninthesystem. Thecorrelationof
higher temperature below T we observe large fluctua- neuronsdecaysmuchfasterinshortrangedinteractionas
c
tions inthe evolutionpatterns withverysmallglobalor- comparedtolong ranged,but itscaleswith r/Lshowing
dering; instead of picking one of the up-spin, down-spin, the universality of neuron interaction in brain.
6
Given the current focus on the biological network and Lett. 29, 917 (1972).
theirfunctionality,wehopethatthispaperwillmotivate [26] N.G.deBruijn,AsymptoticMethodsinAnalysis(Dover,
fresh interestin the evolutiondynamics andmorphology New York,1981).
[27] K. Binder and D. W. Heermann, Monte Carlo Simula-
of active and passive neurons. These kinetic processes
tions in Statistical Physics: An Introduction, Springer-
play an important role in determining the functionality
Verlag, Berlin (1988).
of brain. We emphasize that one can gain a good un-
[28] M. E. J. Newman and G. T. Barkema, Monte Carlo
derstandingoftherelevantneurondynamics(wiringand Methods in Statistical Physics, Oxford University Press,
rewiring inside the brain network) from simple coarse- Oxford (1999).
grained models of the type discussed here. Convincing [29] A.Singh,S.Ahmad,S.Puri,andS.Singh,Eur.Phys.J.
to the fact that re-learning help us to memorize things E, 37, 14002, (2013)
[30] Y. Oono and S. Puri, Phys.Rev. Lett.58, 836 (1987).
forlongerduration. Thusitisanattempttopredictand
[31] Y. Oono and S. Puri, Phys.Rev. A 38, 434 (1988).
an outlook to understand the functionaliies of the brain.
[32] S. Puri and Y. Oono, Phys.Rev. A 38, 1542 (1988).
[33] T. Ohta, D. Jasnow, K. Kawasaki, Phys. Rev. Lett. 49,
Acknowledgments 1223 (1982).
[34] A. J. Bray, S. Puri, Phys.Rev. Lett.67, 2670 (1991).
JGacknowledgesthe fellowshipandRKBSfor providing
[35] G. Porod, in Small-Angle X-Ray Scattering, O. Glatter
financial support from the CSIR, India, under sanction
andO.Kratky(eds.),(AcademicPress,NewYork,1982).
no. 25(0221)/13/EMR-II. [36] Y. Oono and S. Puri, Mod. Phys. Lett. B 2, 861 (1988).
[37] V.M.Eguiluz, D.R.Chialvo, G.A.Cecchi, M.Balikiand
A. Vania. Phys.Rev. Lett.94, 018102 (2005).
[38] A. J. Bray, Phys. Rev.E 47, 3191 (1993).
[39] M. Fu,X. Yu,J. Lu and Y. Zuo. Nature483, 92 (2012).
∗ [email protected] [40] P. Monceau, M. Perreau, and F. Hebert, Phys. Rev. B
[1] Y. Liu, and K.A. Dahmen, Phys. Rev. E 79, 061124 58, 6386 (1998).
(2009).
[2] C. C. Lo, R. P. Bartsch, and P. C. Ivanov, Eur. Phys.
Lett.102, 10008 (2013).
[3] D.R. Chialvo, Physica A 340, 756-765 (2004).
[4] S.PuriandV.Wadhawan(eds.),KineticsofPhaseTran-
sitions (CRC Press, Boca Raton, 2009).
[5] A.J. Bray,Adv.Phys. 43, 357 (1994).
[6] A. Onuki, Phase Transition Dynamics (Cambridge Uni-
versity Press, Cambridge, 2002).
[7] R.J. Glauber, J. Math. Phys., 4, 294 (1963).
[8] I.M. Lifshitz, Zh.Eksp. Teor. Fiz. 42, 1354 (1962).
[9] S. M. Allen and J. W. Cahn, Acta Metall. 27, 1085
(1979).
[10] K.Kawasaki, Phys.Rev., 145, 224 (1966).
[11] A.J. Bray,Phys. Rev.Lett. 62, 2841 (1989).
[12] A.J. Bray,Phys. Rev.B41, 6724 (1990).
[13] J.M.Beggs,andD.Plenz,J.Neurosci.23,11167(2003).
[14] J. E. S. Socoloar, and S. A. Kauffman, Phys. Rev. Lett.
90, 4 (2003).
[15] O.Kinouchi, and M. Copelli, Nat. Phys. 2, 348 (2006).
[16] G.Toulouse,S.Dehaene,andJ.P.Changeux,Proc.Natl.
Acad.Sc. 83, 1695 (1986).
[17] D. Fraiman, P. Balenzuela, J. Foss, and D. R. Chialvo,
Phys.Rev.E 79, 061922 (2009).
[18] E. Schneidman, M. J. Berry, R. Segev and W. Bialek,
Nature440, 1007 (2006).
[19] D. Sk. Bassett, A. Meyer-Lindenberg, S. Achard, T.
Duke,andE.Bullmore,Proc.Natl.Acad.Sc.103,19518
(2006).
[20] S. L. Bressler, and V. Menon, Trends in Cognitive Sc.
14, 277 (2010).
[21] A.Schnitzler,andJ.Gross,NatureRev.Neurosc.6,285
(2005).
[22] S.K. Lamoreaux I, Am.J. Phys.67, 850 (1999).
[23] B.J.Hiley,andG.S.Joyce,Proc.Phys.Soc.493(1965).
[24] S.A.Cannas,andF.A.Tamarit,Phys.Rev.B54,12661
(1996).
[25] M. E. Fisher, S.-k. Ma, and B. G. Nickel, Phys. Rev.