An Ising model can be specified as a neural network of N binary variables only constrained by pairwise correlations. Units can have a value of +1 or −1 and are affected by local bias hi and couplings Jij between pairs of units. These parameters take continuous values and we assume couplings to be symmetric with Jii = 0. In order to simulate the behaviour of the model, units are updated sequentially in a random order using Glauber dynamics, by which each unit is activated with a probability that follows a sigmoid function:1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({s}_{i}(t+1))={[1+{e}^{-\beta 2{H}_{i}(t){s}_{i}(t+1)}]}^{-1}$$\end{document}P(si(t+1))=[1+e−β2Hi(t)si(t+1)]−1where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${H}_{i}(s)={h}_{i}+{\sum }_{j}{J}_{ij}{s}_{j}(t)$$\end{document}Hi(s)=hi+∑jJijsj(t) is the effective field received by neuron i summing the internal field and inputs from other neurons (and 2Hi(s)si(t + 1) is the energy difference required to flip the sign of unit i). The state of the model will be updated by sequentially applying Glauber dynamics (i.e. Equation 1) each simulation step to all units of the network in a random order. When updated sequentially, an Ising model with symmetric couplings will reach an equilibrium maximum entropy distribution:2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(s)=\frac{1}{Z}{\exp }[\beta (\sum _{i}{h}_{i}{s}_{i}+\sum _{i < j}{J}_{ij}{s}_{i}{s}_{j})]$$\end{document}P(s)=1Zexp[β(∑ihisi+∑i<jJijsisj)]where the distribution follows an exponential family \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(s)=\frac{1}{Z}{e}^{-\beta E(s)}$$\end{document}P(s)=1Ze−βE(s) and Z is a normalization value. The energy E(s) of each state of size N is defined in terms of the bias hi and couplings Jij between pairs of units, with β = 1/(TkB), being kB the Boltzmann’s constant and T the temperature of the system. Without loss of generality we can set an operating working temperature such that β = 1.