# The effect of a density dependent bag constant on the structure of hot
neutron star with a quark core

T. Yazdizadeh1 and G. H. Bordbar2,3 1Islamic Azad University, Bafgh Branch,
Bafgh, Iran
2Department of Physics, Shiraz University, Shiraz 71454, Iran (e-mail:
<EMAIL_ADDRESS>
3Research Institute for Astronomy and Astrophysics of Maragha,
P.O. Box 55134-441, Maragha 55177-36698, Iran

###### Abstract

As we go from center toward the surface of a neutron star, the state of
baryonic matter changes from the de-confined quark-gluon to a mixed phase of
quark and hadronic matter, and a thin crust of hadronic matter. For the quark
matter, within MIT bag model, the total energy density of the system is the
kinetic energy for non-interacting quarks plus a bag constant. In this article
first we have considered a density dependent bag constant obtained using the
recent experimental results of CERN SPS on the formation of a quark-gluon
plasma. For calculations of the hadron phase, we use the lowest order
constrained variational method. The equation of state of mixed phase has been
determined using Gibbs conditions. Finally, we have calculated the structure
of a hot neutron star with quark core employing TOV equation. Our results show
that a density dependent bag constant leads to a higher mass and lower radius
for the hot neutron star with respect to the case in which we use a fixed bag
constant.

##  I Introduction

Neutron stars (NS) are among the densest of massive objects in the universe. A
hot neutron star is born following the gravitational collapse of the core of a
massive star just after the supernova explosion. The interior temperature of a
neutron star at its birth is of the order $20-50\ MeV$ burrows1986 . As we go
from the surface to the center of a neutron star, at sufficiently high
densities, the matter is expected to undergo a transition from hadronic matter
where the quarks are confined inside the hadrons to a state of deconfined
quarks, with up, down and strange quarks in the quark matter. Other quarks
have higher masses and do not appear in this state. Glendenning has shown that
a proper construction of the hadron-quark phase transition inside the neutron
stars implies the coexistence of nucleonic matter and quark matter over a
finite range of pressure. Therefore, a mixed hadron-quark phase exists in the
neutron star where its energy is lower than that of the quark matter and the
nucleonic matter over a finite rang of Pressure glendenning1992 . This shows
that we can consider a neutron star to be composed of a hadronic matter layer,
a mixed phase of quarks and hadrons and, in the core, quark matter.

The deconfined quark phase is treated within the popular MIT bag model
chodos1974 . In this model, the total energy density is the sum of a
nonperturbative energy shift ${\cal B}$ (the bag constant) and the kinetic
energy for noninteracting quarks. The bag constant ${\cal B}$ can be
interpreted as the difference between the energy densities of the perturbative
vacuum and physical ones, which has a constant value such as ${\cal B}=55$ and
$90MeV/fm^{3}$ in the initial model of MIT, recently it is constrained to be
compatible with the recent experimental results obtained at CERN on the
formation of a quark-gluon plasma. The resulting picture is the following:
during the early stages of the heavy-ion collision, a very hot and dense state
(fireball) is formed whose energy materializes in the form of quarks and
gluons strongly interacting with each other, exhibiting features consistent
with expectations from a plasma of deconfined quarks and gluons rafelsky1982 ;
matsui1986 . In general, it is not obvious if the informations on the nuclear
EOS from high energy heavy ion collisions can be related to the physics of
neutron stars interior. The possible quark-gluon plasma produced in heavy ion
collision is expected to be characterized by small baryon density and high
temperature, while the possible quark phase in neutron stars appears at high
baryon density and low temperature. However, if one adopts for the hadronic
phase a non-interacting gas model of nucleons, antinucleons and pions, the
original MIT bag model predicts that the deconfined phase occurs at an almost
constant value of the quark-gluon energy density, irrespective of the
thermodynamical conditions of the system cleman1986 .

Burgio et al. have investigated the structure of neutron stars with a quark
core at zero Burgio2002a  and finite temperatures burgio2007 , using the
Brueckner-Bethe-Goldstone formalism to determine the equation of state of the
hadronic matter. We have calculated the structure properties of the cold
neutron star by considering a quark phase at its core bordbar2006  and
compared the results with our previous calculations for the neutron star
without the quark core bordbarh2006 . In these works, we have employed the
lowest order constrained variational (LOCV) method for the hadronic matter
calculations. Recently, we have calculated the structure of a hot neutron star
with a quark core with a fixed bag constant (${\cal B}=90\ MeV/fm^{3}$)
yazdizadeh2011 . In the present paper, we intend to extend these calculations
to hot neutron stars with a quark cores by considering a density dependent bag
constant.

##  II  Equation of State

In this section, we calculate the equation of state of a neutron star composed
of a hadronic matter, a mixed phase of quarks and hadrons and a quark core as
follows.

###  II.1  Hadron Phase

In our calculations, the equation of state of hot nucleonic matter is computed
using the lowest order constrained variational (LOCV) method bordbar1997 ;
bordbar1998 ; modarres1998 ; bordbar2007a ; bordbar2007b ; bordbar2008a ;
bordbar2008b ; bigdeli2009 . The details of our calculations for the hadronic
phase has been fully discussed in yazdizadeh2011 .

###  II.2 Quark Phase

We use the MIT bag model for the quark matter. In this model, the energy
density is the sum of kinetic energy of quarks and a bag constant (${\cal B}$)
which is interpreted as the difference between energy densities of non
interacting quarks and interacting ones farhi1984 ,

${\cal E}_{tot}={\cal E}_{u}+{\cal E}_{d}+{\cal E}_{s}+{\cal B},$ (1)

where ${\cal E}_{i}$ is the kinetic energy per volume of particle $i$, the
kinetic energy of quarks has been discussed in yazdizadeh2011  and The bag
constant ${\cal B}$, can be interpreted as the difference between the energy
densities of the noninteracting quarks and interacting ones, which has
constant values such as ${\cal B}=55,\ 90$ and $220\ MeV$ in the initial model
of MIT. But the density of quark matter is not uniform, therefor we can
consider a density dependent ${\cal B}$. We try to determine a range of
possible values for ${\cal B}$ by exploiting the experimental data obtained at
the CERN SPS, where several experiments using high-energy beams of Pb nuclei
reported (indirect) evidence for the formation of a quark-gluon plasma
jacob2000 ; heinz2000 . According to the analysis of those experiments, the
quark-hadron transition takes place at about seven times normal nuclear matter
energy density $({\cal E}_{0}=156MeVfm^{-3})$. In the literature, there are
attempts to understand the density dependence of ${\cal B}$ adami1993 ;
blaschke1999 . However, currently the results are highly model dependent and
no definite picture has come out yet. Therefore, we attempt to provide
effective parametrization for this density dependence. Our parametrization are
constructed in such a way that at asymptotic densities ${\cal B}$ approaches a
finite value ${\cal B}_{\infty}$. For the bag constant (${\cal B}$), we use a
density dependent Gaussian parametrization Burgio2002a ; Baldo2006 ,

${\cal B}(n)={\cal B}_{\infty}+\left({\cal B}_{0}-{\cal
B}_{\infty}\right)\exp[-\beta(n/n_{0})^{2}]$ (2)

with ${\cal B}_{0}={\cal B}(n=0)=400\ MeV/fm^{3}$ and $\beta=0.17$. We know
that the value of the bag constant (${\cal B}$) should be compatible with
empirical results. The experimental results at CERN-SPS show a proton fraction
$x_{p}=0.4$ Burgio2002b . Therefore, for calculation of ${\cal B}_{\infty}$,
we employ the equation of state of the asymmetric nuclear matter as follows.
First, we use the equation of state of asymmetric hadronic matter
characterized by a proton fraction $x_{p}=0.4$ and the UV14 + TNI potential.
By assuming that the hadron-quark transition takes place at the energy density
${\cal E}=1100MeVfm^{-3}$, we find that the hadronic matter baryonic density
is $n_{t}$ (transition density) and at values lower than $n_{t}$, the quark
matter energy density is higher than that of nuclear matter, while by
increasing the baryonic density, the two energy densities become equal at this
density, and after that, the nuclear matter energy density remains always
higher. Energy density equation of quark matter with two flavors u and d
reduces to

${\cal E}_{Q}={\cal E}_{u}+{\cal E}_{d}+{\cal B}$ (3)

Beta-equilibrium and charge neutrality conditions lead to the following
relation for the number density of quarks,

$n_{B}=2n_{u}=n_{d}.$ (4)

We determine ${\cal B}_{\infty}$ by putting quark energy density and hadronic
energy density equal to each other at any temperature. Later we can calculate
the energy density and determine pressure and the equation of state for quark
phase such as determined in yazdizadeh2011 .

###  II.3  Mixed phase

For mixed phase, where the fraction of space occupied by quark matter smoothly
increase from zero to unity, we have a mixture of hadrons, quarks and
electrons. According to the Gibss equilibrium condition, the temperatures,
pressures and chemical potentials of both phases are equal glendenning1992 .
Calculation of equation of state of mixed phase has been fully discussed in
yazdizadeh2011 .

###  II.4 Results

Our results for the energy densities corresponding to different phases i.e.
hadronic matter, pure quark matter and a mixed phase of quarks and hadrons,
are given in Fig. 1 at two temperatures. It can be seen that at low densities
there is a pure hadronic phase. We have found that a mixed phase exists at
higher densities up to $n\sim 0.5fm^{-3}$. It is obvious that a pure quark
phase emerges by increasing the density. The energy density for neutron star
with density dependent and density independent ${\cal B}$ are shown in Fig. 2.
It is clear that for a density dependent bag constant, the energy density is
smaller than that of the independent one. We have calculated the equation of
state (the pressure versus baryon density) for the neutron star with the quark
core using the density dependent and independent bag constants. Fig. 3 shows
the relevant results at two temperatures. These equations of state are used as
an input into the general relativistic equation of hydrostatic equilibrium.

##  III Structure of the Hot Neutron Star with a Quark Core

The equilibrium energy density distribution of slowly rotating spherical star
is determined by the Tolman-Oppenheimer-Volkoff equation (TOV) shapiro1983 ;
glendenning2000 ; weber1999 ,

$\frac{dP}{dr}=-\frac{G[{\cal E}(r)+\frac{P(r)}{c^{2}}][m(r)+\frac{4\pi
r^{3}P(r)}{c^{2}}]}{r^{2}[1-\frac{2Gm(r)}{rc^{2}}]},$ (5)

$\frac{dm}{dr}=4\pi r^{2}{\cal E}(r).$ (6)

$P$ is the pressure and ${\cal E}$ is the total energy density. For a given
equation of state in the form $P({\cal E})$, the TOV equation yields the mass
and radius of star as a function of central energy density.

In our calculation for the hot neutron star with quark core, we use the
following equation of state: (i) Below the density $0.05fm^{3}$, we use the
equation of state calculated by Baym baym1971 . (ii) For the hadron phase,
from the density of $0.05fm^{3}$ up to the density where the mixed phase is
started, we use the equation of state which is calculated in section II.1.
(iii) In the range of densities in mixed phase, we use the equation of state
calculated in section II.3. (iv) For quark phase, we use the equation of state
calculated in section II.2. Using the above equations of state, we integrate
the TOV equation numerically and determine the structure of this star. All
calculations are done for the density dependent bag constant (${\cal B}(n)$)
at two different temperatures $T=10$ and $20MeV$. Our results are as follows.

The gravitational mass versus the central mass density for a hot neutron star
with quark core for two different temperatures has been presented in Fig. 4.
We can see there is the limiting mass for hot neutron star and this mass
increases when we consider a density dependent bag constant in MIT model. This
is reasonable because when a density dependent bag constant is considered, the
equation of state is softer as this is seen in Fig. 3 for two temperatures.
The radius as a function of central mass density for hot neutron star with
quark core for two different temperatures has been presented in Fig. 5. The
radius of star decreases when the bag constant is density dependent. Our
results for the maximum gravitational mass of the hot neutron star with the
quark core and the corresponding values of radius and central mass density
have been given in Tables 1 and 2 for two different temperatures.

##  IV Summary and Conclusion

The structure of hot neutron star with a quark core using a density dependent
bag constant has been investigated. From the surface toward the center of hot
neutron star, a pure hadronic matter, a mixed phase of quarks and hadrons in a
range of densities determined by employing the Gibbs conditions, and a pure
quark matter in the core, have been considered. We have employed the LOCV
method at finite temperature to get the equation of state of hot hadronic
matter. The MIT bag model with the density dependent bag constant obtained by
using the recent experimental results of CERN SPS on the formation of a quark-
gluon plasma has been applied to compute the equation of state of hot quark
matter. We have solved the TOV equation by a numerical method to determine the
structural properties of hot neutron star with the quark core at T = 10 and 20
MeV. The results have been compared to those for the hot neutron star with
${\cal B}=90\ MeV/fm^{3}$. We have found that a density dependent bag constant
leads to a higher mass and a lower radius for the neutron star in comparison
to the case in which the constant ${\cal B}=90\ MeV/fm^{3}$ has been used.

## Acknowledgements

We wish to tank the Research Council of Islamic Azad University, Bafgh Branch.
G. H. Bordbar wishes to thank the Research Institute for Astronomy and
Astrophysics of Maragha and Shiraz University Research Council for financial
support.

## References

  * (1) A. Burrows and J. M.Lattimer,_ApJ_ 307, 178 (1986)
  * (2) N. K. Glendenning, _Phys. Rev._ D46, 1274 (1992)
  * (3) A. Chodos, R.L. Jaffe, K. Johnson, C. B. Thorn and V. F. Weisskopf, _Phys. Rev._ D9, 3471( 1974).
  * (4) J. Rafelsky and B. Muller, _Phys. Rev. Lett._ 48 1066(1982)
  * (5) T. Matsuiand and H.Satz, _Phys. Lett._ B178, 416(1986)
  * (6) J. Cleymans, R.V. Gavai and E. Suhonen, _Physics Rep._ 130, 217 (1986).
  * (7) G. F. Burgio, M. Baldo, P. K. Sahu and H.-J. Schulze, _Phys. Rev._ C66, 025802(2002)
  * (8) G. F. Burgio, M. Baldo, O. E. Nicotra and H.-J. Schulze,_Ap &SS _308, 387(2007)
  * (9) G. H. Bordbar, M. Bigdeli and T. Yazdizadeh, _Int. J. Mod. Phys._ A21, 5991 (2006)
  * (10) G. H. Bordbar and M. Hayati, _Int. J. Mod. Phys._ A21, 1555 (2006)
  * (11) T. Yazdizadeh and G.H. Bordbar, _Res. Astron. Astrophys_ 11, 471(2011)
  * (12) G. H. Bordbar and M. Modarres, _J. Phys. G: Nucl. Phys._ 23, 1631(1997)
  * (13) G. H. Bordbar and M. Modarres, _Phys. Rev._ C57, 714.(1998)
  * (14) M. Modarres and G. H. Bordbar, _Phys. Rev._ C58, 2781(1998)
  * (15) G. H. Bordbar and M. Bigdeli, _Phys. Rev._ C75, 045804(2007)
  * (16) G. H. Bordbar and M. Bigdeli, _Phys. Rev._ C76, 035803(2007)
  * (17) G. H.Bordbar and Bigdeli _Phys. Rev._ C77, 015805(2008)
  * (18) G. H.Bordbar and Bigdeli _Phys. Rev._ C78, 054315(2008)
  * (19) M. Bigdeli, G. H. Bordbar and Z. Rezaei,_Phys. Rev._ C80, 034310 (2009).
  * (20) E.Farhi, R. L. Jaffe, _Phys. Rev._ D30, 2379(1984)
  * (21) U.Heinz and M. Jacobs, nucl-th/0002042.
  * (22) U.Heinz, _Nuch. phys._ A685, 414(2001)
  * (23) C. Adami and G.E. Brown, _Phys. Rep._ 234, 1(1993)
  * (24) D. Blaschke, H. Grigorian, G. Poghosyan, C.D. Roberts and S. Schmidt, _Phys. Lett._ B450, 207(1999)
  * (25) M. Baldo, G. F. Burgio and H.-J. Schulze, _SUPERDENSE QCD MATTER AND COMPACT STARS NATO Science Series_ Volume 197, Part II, 113 (2006)
  * (26) G. F. Burgio, M. Baldo ,P. K. Sahu, A. B. Santra and H.-J. Schulze, _Phys. Lett._ B526, 19 (2002)
  * (27) S. L. Shapiro and S. A. Teukolski, 1983, _Black Holes, White Dwarfs and neutron Stars_ (New York ).
  * (28) N. K. Glendenning, _Compact Star, Nuclear Physics, Particle Physics, and General Relativity_ (Springer, New York), 2nd ed.(2000)
  * (29) F. Weber, _Pulsars as Astrophysical Laboratories for Nuclear and Particle Physics_ (Institute of Physics, Bristol)(1999)
  * (30) G. Baym, C. Pethick and P. Sutherland, _Astrophys. J._ 170, 299(1971)

Table 1: Maximum gravitational mass $(M_{max})$ and the corresponding radius $(R)$ and central mass density $({\cal E}_{c})$ of hot neutron star with a quark core for density dependent and fixed ${\cal B}$ at $T=10\ MeV$. | $NS+QuarkCore$ | $M_{max}\left(M_{\,\odot}\right)$ | $R\left(km\right)$ | ${\cal E}_{c}\left({10^{14}gr/cm^{3}}\right)$
---|---|---|---|---
| Dependent ${\cal B}$ | 2.032 | 10.39 | 24.42
| Fixed ${\cal B}$ | 1.76 | 10.45 | 27.38
Table 2: Maximum gravitational mass $(M_{max})$ and the corresponding radius $(R)$ and central mass density $({\cal E}_{c})$ of hot neutron star with a quark core for density dependent and fixed ${\cal B}$ at $T=20\ MeV$. | $NS+QuarkCore$ | $M_{max}\left(M_{\,\odot}\right)$ | $R\left(km\right)$ | ${\cal E}_{c}\left({10^{14}gr/cm^{3}}\right)$
---|---|---|---|---
| Dependent ${\cal B}$ | 2.033 | 10.9 | 24.43
| Fixed B | 1.78 | 11 | 27.37

Figure 1: Energy density versus the baryon density for the hadron phase (solid
line), mixed phase (dashed line) and quark phase (dotted line) at $T=10$ (a)
and $20\ MeV$ (b).

Figure 2: Energy density versus the baryon density for the neutron star with
the quark core with density dependent (dotted line) and independent (solid
line) bag constant at $T=10$ (a) and $20\ MeV$ (b).

Figure 3: Pressure versus the baryon density for the neutron star with the
quark core with density dependent (dotted line) and independent (solid line)
bag constant at $T=10$ (a) and $20\ MeV$ (b).

Figure 4: Gravitational mass versus the central mass density for the neutron
star with the quark core with density dependent (dotted line) and independent
(solid line) bag constant at $T=10$ (a) and $20\ MeV$ (b).

Figure 5: radius versus the central mass density for the neutron star with the
quark core with density dependent (dotted line) and independent (solid line)
bag constant at $T=10$ (a) and $20\ MeV$ (b).


Evaluation of dynamic causal modelling and Bayesian model selection using
simulations of networks of spiking neurons

Matthew G. Thomas

Anaesthetic Department, Bedford Hospital, United Kingdom

Former affiliation: School of Medicine and Biomedical Sciences, University of
Oxford, United Kingdom

<EMAIL_ADDRESS>
Inferring the mechanisms underlying physiological and pathological processes
in the brain from recorded electrical activity is challenging. Bayesian model
selection and dynamic causal modelling aim to identify likely biophysical
models to explain data and to fit the model parameters. Here, we use data
generated by simulations to investigate the effectiveness of Bayesian model
selection and dynamic causal modelling when applied at steady state in the
frequency domain to identify and fit Jansen-Rit models. We first investigate
the impact of the necessary assumption of linearity on the dynamics of the
Jansen-Rit model. We then apply dynamic causal modelling and Bayesian model
selection to data generated from simulations of linear neural mass models,
non-linear neural mass models, and networks of discrete spiking neurons.
Action potentials are a characteristic feature of neuronal dynamics but have
not previously been explicitly included in simulations used to test Bayesian
model selection or dynamic causal modelling. We find that the assumption of
linearity abolishes the qualitative transitions seen as a function of the
connectivity parameter in the original Jansen-Rit model. As with previous
work, we find that the recovery procedures are effective when applied to data
from linear Jansen-Rit neural mass models, however, when applying them to non-
linear neural mass models and networks of discrete spiking neurons we find
that their effectiveness is significantly reduced, suggesting caution is
required when applying these methods.

Introduction
Electroencephalography (EEG) provides a non-invasive method for spatially
resolved measurement of the electrical activity of the brain1. EEG is a
valuable clinical and research tool, however, it remains challenging to infer
key properties of the underlying network, such as channelopathies, from EEG
recordings.

Dynamic causal modelling2, 3, 4, 5, 6, 7, 8, 9 (DCM) aims to infer the
parameters of biophysical models from recorded electrical activity.
Biophysical models are constructed such that model parameters have a direct
interpretation in terms of network properties, for example, connection weights
or ion channel dynamics. By inferring the parameters of biophysical models
from EEG or local field potential (LFP) recordings, DCM may provide a means to
infer properties of the underlying neuronal network. While DCM seeks to infer
the parameters of a particular biophysical model, choosing a particular model
requires knowledge of the structure of the underlying network. Bayesian model
selection (BMS)3, 10, 8 seeks to infer which of a range of models best
explains the recorded data.

To date, validation of these methods has taken two forms. Firstly, the
techniques have been applied to data generated from models with the same
structure as those underlying the DCM and BMS techniques and the recovered
parameters were compared to those used as input to the original model 2, 3, 5,
11, 12, 13. The second approach to validation has been the application to data
from human subjects and animals2, 4, 5, 11, 14, 10, 15, 16, 17, 18.

The use of simulations in validation has the advantage of exact knowledge of
the network structure and parameters producing the synthetic data, allowing
direct comparison with those recovered. To date, however, this approach has
been applied by generating data with equations of exactly the form employed in
the recovery methods2, 3, 5, 11, 12, 13. This is likely to have overestimated
their effectiveness. Furthermore, the model equations that have been used to
generate synthetic data have had no explicit spiking dynamics and thereby omit
a characteristic feature of the dynamics of networks of neurons.

Studies employing data from humans and animals have been of two types. Some
studies used patients with previously diagnosed channelopathies 19, 17. The
strength of this approach is that the underlying pathology is known and thus
can be used to validate the results of the analysis methods. The weakness,
however, is the extremely limited number of patients and channelopathies. The
alternative approach is to study evoked responses or heterogeneous disease
groups. Much larger sample numbers can be obtained, however, the true network
properties are not completely certain, making validation of the results
difficult 2, 4, 5, 11, 14, 10, 15, 16, 17, 18.

Here, we further investigate the validity of DCM and BMS, considering their
implementation based on Jansen-Rit models in frequency space and applied at
steady state5. We start by analysing the effects of the assumptions used to
generate the necessary transfer functions. We then apply DCM and BMS to
simulated data. We move from simulations of linear neural mass models, through
non-linear neural mass models, to simulations of networks of discrete neurons
with spiking dynamics. This work builds on past computational studies, with
their control and exact knowledge of the network structure and parameters to
be recovered, by adding spiking activity, a key feature of neuronal function
not implemented in previous forward models and a key difference between the
underlying system of interest and the neural mass models that DCM and BMS fit
to data.

The use of transfer functions to fit models to data in frequency space
requires an assumption of linearity5. We find that implementing this linearity
in Jansen-Rit models causes a qualitative change in the model dynamics, losing
the transitions otherwise seen as a function of coupling parameters. As in
previous work5, we find good agreement between the models and parameters
recovered from data generated from linearised Jansen-Rit models, however, we
find that the introduction of non-linearity in the neural mass models and
networks of discrete neurons significantly reduces the effectiveness of both
DCM and BMS.

Jansen-Rit Neural Mass Models
Neural mass models simplify the complexity of networks of neurons to more
tractable averaged quantities20, 21. The Jansen-Rit model22, 23 is a model of
a local cortical circuit. It consists of components reflecting three
populations of neurons: one population of feedforward pyramidal neurons and
two populations of interneurons, one excitatory and one inhibitory, that form
local feedback loops (Figure 1A). The feedforward pyramidal neurons receive
external input and give external output.

Figure 1: Model structures. (A) The Jansen-Rit neural mass model, composed of
populations of feedforward pyramidal cells, excitatory interneurons, and
inhibitory interneurons. $P(t)$ is the external input; $h_{e}(t)$ and
$h_{i}(t)$ are the excitatory and inhibitory impulse response functions
respectively; $Sigm$ are sigmoid functions; coupling constants are as defined
in the text; and addition and subtraction symbols indicate excitation and
inhibition respectively. (B) Network of spiking neurons. Connections reflect
the structure of the Jansen-Rit model. Black lines represent excitatory
connections, red lines represent inhibitory connections. The directions of the
connections are indicated by arrows. A reduced number of neurons in each
population is shown for clarity.

Each population in the model can be thought of as being composed of two
blocks. The first block takes a mean input firing rate and outputs a mean
membrane potential. The impulse response of this block has the form of a post-
synaptic potential, for $t\geq 0$: $aAte^{-at}$ for an excitatory input spike
and: $bBte^{-bt}$ for an inhibitory input spike, where $A$ and $B$ are
amplitude parameters and $a$ and $b$ are reciprocals of the time constants of
the post-synaptic potential. The second block takes the mean membrane
potential from the first block and produces a mean output spike rate. This
transformation takes the form of a sigmoid:

$S(v)=\frac{v_{Max}}{1+e^{r(v_{0}-v)}}-\frac{v_{Max}}{1+e^{rv_{0}}}$ (1)

where $r$ controls the gradient of the sigmoid function, $v_{Max}$ reflects
the maximum firing rate of the population, and $v_{0}$ reflects the potential
at which the firing rate is half of maximum. The second term in the expression
shifts the function so that $S(0)=0$.

The external input, $P(t)$, is white noise. Here it is uniformly distributed
across the range from 120 to 320 pulses per second. The parameters $C_{1}$,
$C_{2}$, $C_{3}$, and $C_{4}$ control the strength of coupling between the
populations and the variables $y_{0,a}$, $y_{0,b}$, $y_{1}$, and $y_{2}$
correspond to the average membrane potentials within the blocks, as in Figure
1A. The output from the model is the average potential of the feedforward
pyramidal neurons: $y_{1}-y_{2}$. Illustrative parameter values are presented
in Table 1.

Parameter | Example value
---|---
$a$ |  $100$ s-1
$A$ | 3.25 mV
$b$ |  $50$ s-1
$B$ | 22 mV
$C_{1}$ | 300
$C_{2}$ | $0.8C_{1}$
$C_{3}$ | $0.25C_{1}$
$C_{4}$ | $0.25C_{1}$
$v_{Max}$ |  $5$ s-1
$v_{0}$ |  $6$ mV
$r$ |  $0.56$ mV-1

Table 1: Illustrative Jansen-Rit model parameters.

The model can be described by four coupled second order ordinary differential
equations:

$\frac{d^{2}y_{0,a}(t)}{dt^{2}}=aAC_{1}S(y_{1}(t)-y_{2}(t))-2a\frac{dy_{0,a}(t)}{dt}-a^{2}y_{0,a}(t)$
(2)
$\frac{d^{2}y_{0,b}(t)}{dt^{2}}=aAC_{3}S(y_{1}(t)-y_{2}(t))-2a\frac{dy_{0,b}(t)}{dt}-a^{2}y_{0,b}(t)$
(3)
$\frac{d^{2}y_{1}(t)}{dt^{2}}=aA(P(t)+C_{2}S(y_{0,a}(t))-2a\frac{dy_{1}(t)}{dt}-a^{2}y_{1}(t)$
(4)
$\frac{d^{2}y_{2}(t)}{dt^{2}}=bBC_{4}S(y_{0,b})-2b\frac{dy_{2}(t)}{dt}-b^{2}y_{2}(t)$
(5)

The dynamics of the Jansen-Rit model have been well studied22, 23, 24, 25. The
model produces a range of dynamics depending on the connectivity parameter,
$C$, where $C_{1}=C$. As the connectivity parameter is increased there is a
transition from noisy output to periodic oscillations (Figure 2A).

Figure 2: Dynamics of the Jansen-Rit model. (A) Dynamics with the original
sigmoid relationship between the average membrane potential and output firing
rate. (B) Dynamics with the sigmoid replaced with a linearised relationship.
The relevant coupling constant is shown at the top of each column. (i) Time
courses of $y_{1}-y_{2}$ at steady state. (ii) Phase portraits:
$\dot{y}_{1}-\dot{y}_{2}$ plotted against $y_{1}-y_{2}$. (iii) Amplitudes of
the Fourier transforms of the time courses shown in (i) (black) and the
average of 20 time courses (red). (C) Stability analysis. The poles of the
transfer functions of the linearised models (B) are shown in dark blue. Poles
with positive real values are associated with unstable dynamics.

Dynamic Causal Modelling
DCM aims to infer probability distributions for the parameters of biophysical
models using macroscopic observations2 (Figure 3A); thereby allowing inference
of the mechanistic basis of observed phenomena at the microscopic level from
non-invasive measurements2, 5. DCM has been implemented in a number of forms,
with varying models and recording modalities 26, 2, 3, 4, 11, 5, 27, 6, 10, 7,
18, 9, 8.

Here we will consider the implementation used in Moran et al5 as it is
conceptually simple and so readily allows identification of its strengths and
limitations. This approach uses a Jansen-Rit model as the biophysical model
whose parameters are to be inferred. The data considered are
electrophysiological recordings, such as EEG or LFP time series. Probability
distributions for the parameters of the model are inferred from the data in
frequency space. The transformation to frequency space for the observed data
is readily obtained by taking the Fourier transform of the time series. The
analytic transformation to the form for the model in frequency space requires
the equations of the model to be linear. To achieve this, the sigmoid
functions in the Jansen-Rit model are approximated by linear functions, with
proportionality constant $\gamma$:

$S(v)=\frac{v_{Max}re^{rv_{0}}}{(1+e^{rv_{0}})^{2}}v=\gamma v$ (6)

With this assumption, the transfer function of the model can be found by
taking the Laplace transform of the equations. The transfer function for the
Jansen-Rit model equations in the previous section is readily derived:

$H(s)=\left(\frac{(s+a)^{2}}{aA}-\frac{aAC_{1}C_{2}\gamma^{2}}{(s+a)^{2}}+\frac{bBC_{3}C_{4}\gamma^{2}}{(s+b)^{2}}\right)^{-1}$
(7)

In frequency space the model can be thought of as a filter operating on the
spectra of the input, to produce the spectra of the measured output. We can
write this as: $Y(s)=H(s)U(s)$, where $Y(s)$ and $U(s)$ are the spectra of the
model output and input respectively and $H(s)$ is the transfer function for
the system.

We wish to infer the parameters of the transfer function, and, therefore, the
Jansen-Rit model, from the Fourier transform of data. There are a number of
ways to achieve this. A Bayesian approach can be used to include prior beliefs
about parameter values and variational approaches, such as an expectation-
maximization scheme, can be used to invert the model5.

In our study we wish to treat the simplest cases to identify the effects of
key assumptions on the effectiveness of the DCM approach. To this end we will
consider white noise input to the model, yielding, on average over many
instances: $\langle|Y(s)|\rangle\propto\langle|H(s)|\rangle$. We will also
assume no prior knowledge of the parameters and calculate their likelihood
functions by exhaustive enumeration over a wide range around the true model
parameters.

Assuming information about the phase of the input is unknown, we consider only
the amplitude of the input and output spectra. Taking the amplitude amounts to
combining two components with independent, normally distributed, noise in
quadrature. The result is that the output spectra is the transfer function
multiplied by noise with a Rayleigh distribution.

The likelihood of a given model parameter set, $\theta$, is calculated given
the spectra, $W(f)$, of an experimental recording as follows. We write $x(f)$
as the ratio:

$x(f)=\frac{|W(f)|}{|H(f,\theta)|}$ (8)

We consider the probability density function of the Rayleigh distribution,
with scale parameter $\sigma$:

$R(x,\sigma)=\frac{x}{\sigma^{2}}e^{\frac{-x^{2}}{2\sigma^{2}}}$ (9)

We can then calculate the likelihood of a given parameter set given observed
data with samples at discrete frequencies:

$L(\theta|W)=\prod_{i}R(x(f_{i}))$ (10)

We note that the Rayleigh distribution has a single parameter, $\sigma$, and
that the mean of the distribution can be expressed:
$\sigma\sqrt{\frac{\pi}{2}}$. For data generated by the Jansen-Rit model, if
we allow ourselves knowledge of the amplitude of the input noise we can
utilise the formula for the standard deviation of uniformly distributed white
noise ($\frac{1}{\sqrt{12}}(P_{Max}-P_{Min})$) to divide $W(f)$ by this
amplitude. The likelihoods can then be calculated with
$\sigma=\frac{1}{\sqrt{2}}$ in the Rayleigh distribution. If we allow
ourselves to assume the input noise is white, but do not allow knowledge of
the amplitude, we calculate the likelihood by normalising both the observed
spectra and the transfer function by dividing by their means. In this case the
likelihoods are calculated with $\sigma=\sqrt{\frac{2}{\pi}}$. A similar
approach can be employed for models of discrete spiking neurons, where
treating the input spikes as a Poisson distributed train of delta functions
allows calculation of the scale parameter for the resulting Rayleigh
distributed noise in frequency space, if the input frequency is known, using
the formula $\langle U(s)\rangle\propto\frac{\sqrt{(\pi\langle
P\rangle)}}{2}$. In each case an appropriate correction is required to account
for the convention used with regard to the number of samples in the numerical
Fourier transform of the model output.

Figure 3: Illustration of DCM and BMS. (A) DCM seeks to identify the most
likely set of parameters to explain recorded data. The black circles show
illustrative data points. The three coloured curves show the same function
with different parameter sets. The blue curve corresponds to the parameter set
most likely to explain the data. (B) BMS seeks to identify the model that best
explains measured data sets. The coloured points show illustrative data sets.
The coloured curves show different models. (i) Model 1 best explains the data
sets. (ii) Model 2 best explaines the data sets.

Bayesian Model Selection
BMS aims to infer the most likely model to explain a set of observations
(Figure 3 B). For example, we can imagine administering a drug that has some
effect on the brain, while not knowing the mechanism of action of the drug. We
consider different possibilities, for example it may prolong the duration of
excitatory postsynaptic potentials (EPSPs) or it may reduce the amplitude of
inhibitory postsynaptic potentials (IPSPs) (for simplicity we will imagine a
drug that has a single effect). We imagine having EEG recordings associated
with different doses of drug. Bayesian model selection then infers the most
likely model to explain the data. Given models $M_{i}$ and data $D$, Bayes
theorem gives:

$P(M_{i}|D)=\frac{P(D|M_{i})P(M_{i})}{P(D)}$ (11)

Here we will consider the simple case of constant priors, yielding as a result
the maximum likelihood estimate (MLE):

$P(M_{i}|D)\propto P(D|M_{i})$ (12)

For each model we integrate over all possible parameters sets, $\theta$:

$P(D|M_{i})=\int P(D|M_{i},\theta)P(\theta|M_{i})d\theta$ (13)

Assuming $P(\theta|M_{i})$ is a constant independent of $M_{i}$, this becomes:

$P(M_{i}|D)\propto\int P(D|M_{i},\theta)d\theta$ (14)

Here we will consider the simple case where our set of trial models is
restricted to 11 models, for each of which a single parameter, $d$, acts
linearly on one of the parameters in the transfer function. For example, the
transfer function in the case where the action of the intervention increases
the amplitude of EPSPs, $A$, is:

$M(s,d)=\left(\frac{(s+a)^{2}}{adA}-\frac{adAC_{1}C_{2}\gamma^{2}}{(s+a)^{2}}+\frac{bBC_{3}C_{4}\gamma^{2}}{(s+b)^{2}}\right)^{-1}$
(15)

When calculating the likelihood of each model we treat noise as for DCM. We
can then calculate the relative likelihood of the different models and
identify the one most likely to explain the data set.

Materials and Methods
Simulations and data analysis were carried out using bespoke scripts in
Octave-6.2.0. Simulations were carried out with time steps of 0.1 ms. Systems
were initialised with firing rates and potentials of zero. For each simulation
the time to steady state was established by comparing the Fourier transform of
a shifting 1 s window of the output either to the model transfer function, in
the case where the transfer function was exact for the model, or the final
second of the simulation. The initial period of each simulation affected by
transients was discarded. Simulations were repeated 20 times for each
parameter set.

Results
Effect of linearisation on Jansen-Rit model dynamics

Deriving the transfer function for a Jansen-Rit model requires the sigmoid
functions to be linearised. Working with the transfer function is then
equivalent to inferring the parameters to a modified form of the Jansen-Rit
equations. We start by investigating the effect of this linearisation on the
dynamics of the Jansen-Rit neural mass model.

In the original Jansen-Rit model the sigmoid function bounds the firing rate
between $\frac{-v_{Max}}{1+e^{rv_{0}}}$ and
$v_{Max}\left(1-\frac{1}{1+e^{rv_{0}}}\right)$. Linearisation of the sigmoidal
response removes these bounds on the firing rate. This results in a
significant change in the model dynamics, as shown in Figure 2B. Noisy
dynamics persist for a larger range of values of the connectivity parameter.
The dynamics eventually transition to a growing oscillation. This transition
can be understood in terms of the poles of the transfer function (Figure 2C).
For small enough values of the connectivity parameter all poles of the
transfer function have negative real components, corresponding to decaying
dynamics. As the parameter is increased, poles cross the imaginary axis of the
complex plane so that they have positive real components. This then
corresponds to an unstable system where the output increases with time. The
analytic form for the value of the connectivity parameter at which this
transition occurs for this model is derived in the Appendix. For the parameter
set used here, the transition occurs at approximately $C=606.6$.

Recovery of parameters from a linearised Jansen-Rit model

We wish to use DCM to infer the parameters of neural networks. We start by
looking at the simplest case - applying DCM to data generated from a
linearised Jansen-Rit model, such that the transfer function we use is exactly
that for the model equations used to generate the data. We allow ourselves
exact knowledge of the structure of the network, the form of the dynamics of
the model, and all but one of the model parameters in turn. We then use DCM to
calculate a maximum likelihood estimate for the unknown parameters.

Figure 4: Recovery of the parameters of the linearised Jansen-Rit model using
DCM. Top row: Relative likelihood of parameter values (C = 250). The vertical
blue bars indicate the value of the parameter used in the model. Bottom row:
Relative error in the maximum likelihood estimates of the parameters as a
function of the coupling constant, $C$, used in the model. The black curves
show the results when the input spike rate is known. The dashed red curves
show the likelihoods derived when the input spike rate is unknown.

Twenty time series were generated from the model at each of a range of
connectivity constants, $C$, noting the dependence of the other connection
constants on $C$. The likelihood function was then calculated for each
parameter in turn, while the remaining parameters were set to their true
values. We note that the likelihood function did not have to be calculated for
every parameter as some parameters only appear in the transfer function
(Equation 7) as a product with other particular parameters. For example,
$C_{1}$ and $C_{2}$ only appear in the transfer function as a product and so a
relative change in one cannot be distinguished from a relative change in the
other. In contrast, the constant $A$ appears in the product with $C_{1}$ and
$C_{2}$, but also appears in other terms in the transfer function and so
variation in $A$ can have effects distinguishable from variation in $C_{1}$
and $C_{2}$.

The likelihood functions were calculated for each parameter both with and
without allowing knowledge of the input rate to the model. Figure 4 shows both
the likelihood functions and the relative error in the MLE for illustrative
parameters. The exhaustive results are presented in the Appendix (Figures 9,
10, 13, and 14). The likelihood function for most of the parameters was
unimodal, however, it was bi-modal for $r$, with similar amplitudes at the two
modes. As $C$ was varied, the MLE would usually relate to the true mode, but
would sometimes be the second value (Figure 4C).

For most parameters, the error in the MLE was greater at smaller values of the
connectivity parameter, $C$. This can be understood in terms of the smaller
effect on the dynamics of the system for most parameters when the coupling is
weak. The parameter $a$ could be accurately recovered regardless of the
strength of the connection because it appears in a term of the transfer
function even when $C_{1}=0$. The parameter, $A$, also appears in the transfer
function when $C_{1}=0$ and so can be recovered with weak coupling when the
amplitude of the input to the model is known. However, when the input to the
model is unknown and the likelihood estimate is calculated with normalised
spectra the effect of $A$ is lost and its value cannot be recovered for small
values of the connectivity parameter.

Recovery of parameters from a Jansen-Rit model without linearisation

Neuronal dynamics are inherently non-linear. We are interested in the
effectiveness of DCM at recovering the parameters of such non-linear systems
when the inference procedure inherently employs a linear model.

In moving towards testing DCM on data generated from a spiking neural network
we next apply DCM to time series generated from the original, non-linear,
Jansen-Rit model.

We apply the same approach as in the previous section, but replacing the
linearised sigmoid with the original sigmoid function when generating the
data. As we saw above, exchanging the sigmoid function with a linearised form
has a significant impact on the dynamics of the model. The transfer function
corresponds to the linear function while re-introduction of the sigmoid
function to the Jansen-Rit model generates different dynamics with different
amplitudes and peaks in its frequency spectra (Figure 2A). As the value of
$C_{1}$ is increased, the dynamics of the model deviate away from dynamics
exhibited by the linearised model (Figure 2B). The MLEs resulting from the
application of DCM therefore become significantly less accurate than those
calculated for the linearised model (Figure 5). The exhaustive results are
presented in the Appendix (Figures 11, 12, 13, and 14). For the same reasons
as in the linearised case, when the connectivity constant is small, all but
$a$ and $A$ weakly influence the transfer function. Again, the likelihood
function for $r$ is sometimes multimodal.

Figure 5: Recovery of the parameters of the non-linear Jansen-Rit model using
DCM. Top row: Relative likelihood of parameter values (C = 250). The vertical
blue bars indicate the true value of the parameter used in the model. Bottom
row: Relative error in the maximum likelihood estimates of the parameters as a
function of the coupling constant, $C$, used in the model. The black curves
show the results when the input spike rate is known. The dashed red curve
shows the likelihood derived from the same fit but with both the transfer
function and Fourier transform normalised, reflecting the scenario where the
input spike rate is unknown. Large scale turning points in the error, for
example at $C\sim 140$, correspond to transitions in the dynamics of the
Jansen-Rit model. Smaller scale oscillations in the error, for example the
error in the recovery of $r$ when $C>200$, correspond to shifts between the
two modes of the relative likelihood.

Recovery of parameters from unconnected spiking neurons

Action potentials are a key characteristic of neuronal dynamics. Neural mass
models consider averages over these events, yielding smoothed dynamics. Here
we apply DCM to the effective LFP from a simulated population of unconnected
spiking neurons.

To isolate the effect of spiking dynamics on the effectiveness of DCM we chose
the subthreshold dynamics of our model neurons to match the form of the post-
synaptic potential used in the Jansen-Rit model, with impulse response:
$aAte^{-at}$. The only difference between the dynamics is then the spiking
events. As there is no coupling, we can see from the equation for the sub-
threshold dynamics that there are only two parameters to be recovered: $a$ and
$A$.

We simulated three models: (i) neurons with an infinite threshold, giving
purely subthreshold membrane potential dynamics (Figure 6Ai); (ii) neurons
with a finite threshold potential (5 mV above their resting potential) such
that if it is reached, the neuronal potential returns to its resting potential
linearly over a 4 ms refractory period (Figure 6Aii); (iii) neurons behaving
as in (ii) except that when the threshold is reached the potential follows a
FitzHugh-Nagumo28 form action potential ($\tau=12.5$, $a=0.7$, $b=0.8$) before
returning to its resting potential (Figure 6Aiii).

For each model we simulated 20 instances, each composed of 100 neurons. The
input to each neuron in each instance was a train of Poisson distributed
spikes. For each instance, we generated 100 spike trains, one for each neuron.
For a given instance we used the same set of input trains for each of the
three models. For each instance of each model the effective LFP was generated
by averaging the potential from 100 model neurons. While investigating the
effect of spiking dynamics on the accuracy of DCM we allow ourselves knowledge
of the average input rate to the neurons.

The first model produces an LFP that is identical to the dynamics of a Jansen-
Rit model in which the connection parameters are all set to zero. The transfer
function calculated above fits the spectra exactly (Figure 6B) and the MLE
accurately recovers the model parameters.

Subthreshold neuronal dynamics are known to dominate the dynamics of the LFP1.
The second model is, therefore, the most realistic reflection of the LFP from
spiking neurons. As expected, the resets associated with firing events modify
the frequency spectra, causing them to deviate from the transfer function
(Figure 6B). The error in the MLE of both parameters increases with firing
rate (Figure 6C and D). For the same input rate, the neuronal firing rate
increases as the EPSP amplitude ($A$) is increased and as the reciprocal of
the time constants ($a$) is decreased. As the parameters of the model are
varied, the error in the recovered parameters vary with their impact on the
neuronal firing rates (Figure 6C and D).

The third model is presented for completeness. Including action potentials in
the calculated LFP causes its spectra to deviate substantially from the model
transfer function (Figure 6B). Action potentials may contribute to what is
measured using a local field potential electrode, although the subthreshold
dynamics are thought to dominate1.

Figure 6: Recovery of the parameters from a population of identical uncoupled
spiking neurons using DCM. (A) Illustrative neuronal dynamics without spiking
(i), with resets at spikes (ii), and with FitzHugh–Nagumo model action
potentials (iii). (B) Amplitude of the Fourier transforms of the average
potentials of populations of uncoupled neurons without spiking (blue, 100
neurons in each population, average of 20 populations shown), with resets at
spikes (black), and with action potentials (red). The magenta line shows the
transfer function of the equivalent non-spiking system. The average input rate
to each neuron is 50 Hz (i) and 200 Hz (ii). (C) Error in the recovered EPSP
amplitude for uncoupled neurons with resets at spikes, as a function of
average firing rate (i) and EPSP amplitude used in the model (ii). (D) Error
in the recovered EPSP time constant for uncoupled neurons with resets at
spikes, as a function of average firing rate (i) and EPSP time constant used
in the model (ii). The average input rates to the neurons in Cii and Dii are
100 Hz (black) and 200 Hz (red). The solid black lines indicate the true
parameter values.

Recovery of parameters from networks of spiking neurons

Finally, we investigated the effectiveness of DCM at recovering the parameters
from a simulated network of coupled spiking neurons. We arranged the network
to be as similar as possible to that of the Jansen-Rit model, both in terms of
the structure of the network (Figure 1B) and the inter-spike dynamics of the
neurons. We only included the feedforward population and the inhibitory
interneurons ($C_{1}=C_{2}=0$) as this arrangement results in negative
feedback, yielding a stable simulation.

We simulated networks made up of 750 feedforward pyramidal neurons and 750
inhibitory interneurons. There were no connections within populations. There
were connections between some feedforward neurons and some inhibitory neurons.
Connections were introduced at random. Before modification of parameters, the
connection probability between neurons of both populations was 10% .

The inter-spike dynamics were identical to the Jansen-Rit model response
functions. As in model (ii) in the previous section, neurons had a finite
threshold potential such that if it was reached, the neuronal potential
returned to its resting potential linearly over a 1 ms refractory period 29.

We wished to identify thresholds so that the average firing rates of the
neurons in our model matched those of the populations of the Jansen-Rit model
when the input firing rate to the model matched that of the Jansen-Rit model
at steady state. On this basis, the thresholds of the feedforward neurons and
the inhibitory interneurons were set as 13.52 mV and 4.96 mV respectively. The
external input to the feedforward neurons was a Poisson distributed spike
train with an average rate of 220 Hz.

Networks were simulated for ten seconds. The average potential of the
feedforward neurons was used as an effective LFP1 (Figure 7A). The Fourier
transform of the final five seconds of the effective LFP was used as the
spectrum from each simulation. Networks were simulated 20 times for each
parameter set. The model parameters ($a$, $A$, $b$, $B$, $C_{3}$, and $C_{4}$)
were varied individually from 50% to 200% of the values in Table 1, while the
others were set as in Table 1.

We first compared the spectra of the model LFPs to the transfer functions
generated with the same parameters as used in the models (Figure 7B). The form
of the spectra is similar to that of the transfer function, typically
exhibiting a single maximum, but does not match the shape exactly (Figure 7B
iii). We show illustrative results for the relationship between the spectra
and transfer function as $C_{3}$ is varied in Figure 7C and exhaustive results
in the Appendix (Figure 15). We quantify the relationship using Spearmans rank
correlation coefficient (Table 2). As the parameters are varied the frequency
at which the peak occurs in the spectra is correlated with the frequency at
which the peak occurs in the transfer function for $A$, $b$, and $C_{3}$. The
amplitude of the peak in the spectra is correlated with the amplitude of the
peak in the transfer function for all parameters. The MLE of the parameters is
correlated with the true value of the parameters for $A$ and $C_{3}$.

Parameter | Frequency | Peak | MLE
---|---|---|---
$a$ | -0.45 | 1.00 | -0.32
$A$ | 0.96 | 0.89 | 0.93
$b$ | 0.96 | 0.89 | 0.39
$B$ | -0.66 | 1.00 | -0.18
$C_{3}$ | 0.94 | 1.00 | 0.86
$C_{4}$ | -0.09 | 1.00 | 0.43

Table 2: Spearmans rank correlation coefficients for the relationship between
features of the model spectra and the equivalent transfer function as each
parameter is varied.

Figure 7: Recovery of the parameters of a network of coupled spiking neurons
using DCM. (A) Example membrane potential of a feedforward pyramidal neuron
(i), of an inhibitory neuron (ii), and example LFP from a population of
pyramidal neurons (iii). (B) The transfer function (i) and the spectra of the
LFP of the spiking network (ii) as $C_{3}$, is varied from 50% (red) to 200%
(magenta) of its standard value. (iii) The normalised transfer function (red)
and spectra of the LFP (black) overlaid ($C_{3}=75$). The blue curve shows the
transfer function for the maximum likelihood estimate of $C_{3}$ ($C_{3}=91$).
(C) (i) The relationship between the frequency at the maxima of the transfer
function and the spectra of the LFP as $C_{3}$ is varied. (ii) The
relationship between the amplitude of the maxima of the transfer function and
the spectra of the LFP as $C_{3}$ is varied. (iii) The relationship between
the maximum likelihood estimate of the connectivity from fitting the transfer
function to the spectra and the true value of $C_{3}$. (D) (i) The relative
frequency at the maxima of the transfer function and the spectra of the LFP as
each of the parameters is varied (colours reflect different parameters - full
results in Appendix Figure 15). (ii) The relationship between the relative
amplitude of the maxima of the transfer function and the spectra of the LFP as
each of the parameters is varied. (iii) The relationship between the relative
maximum likelihood estimates of the parameters from fitting the transfer
function to the spectra and the true parameters. In each case the solid red
line shows a least squares fit to the data and the dashed black line shows the
line $y=x$.

Bayesian Model Selection

The effectiveness of BMS was investigated by simulating models and varying
each of the parameters in turn. The model parameters ($a$, $A$, $b$, $B$,
$C_{3}$, and $C_{4}$) were varied individually from 50% to 200% of the values
in Table 1, while the others were set as in Table 1. Each combination of
parameters was simulated 20 times and the Fourier transform calculated as
described above for DCM. BMS was then applied to the data sets for each
parameter.

Application of BMS to data generated from the linearised Jansen-Rit model
correctly identified the model parameter that was varied as one of the most
likely to explain the data in every case. As can be seen from the form of the
transfer function (Equation 7) variation in a number of parameters have
equivalent effects (for example: $B$, $C_{3}$, and $C_{4}$) and so cannot be
distinguished by BMS and it identifies the members of the groups of parameters
as equally likely (Figure 8A).

Applying the same approach to data generated from the Jansen-Rit model without
linearisation and to the model of a network of spiking neurons as described
above found BMS to be less effective (Figure 8B and C), not identifying any of
the parameters correctly in the case of the model of the network of spiking
neurons.

BMS was not applied to the data generated by the model of unconnected spiking
neurons because the model only requires two parameters (EPSP time constant and
amplitude) to describe its dynamics.

Figure 8: Bayesian Model Selection. The rows correspond to the model parameter
that was varied. Within each row the grayscale colour map reflects the
relative likelihood of the simulated data being explained by variation in the
parameter labelling the column. The green circles indicate the elements that
correspond to the true parameter explaining the variation. (A) Data generated
from the linearised Jansen-Rit model. (B) Data generated from the Jansen-Rit
model without linearisation. (C) The effective LFP from the simulation of a
network of spiking neurons.

Discussion
DCM and BMS have the potential to infer the mechanistic basis of key
physiological and pathological processes from recordings of averaged
electrical activity, such as LFP or EEG. We have used simulations to
investigate the effectiveness of DCM and BMS when applied at steady state, in
frequency space, using a transfer function representation of a Jansen-Rit
model. We started by investigating the effect of the linearisation required to
determine a transfer function on the dynamics of the equivalent Jansen-Rit
model. We found that making the assumption of linearity significantly changes
the dynamics of the model. Applying DCM to simulated data from a linearised
model works well when the transfer function is chosen to match the equations
of the linearised model. We still find, however, that the likelihood function
is sometimes multimodal. Some parameters only occur in the transfer function
as a product with particular other parameters and so cannot be resolved
individually. As the coupling between populations is reduced, the accuracy of
the recovery of the parameters of the interneuron populations reduces.
Knowledge of the input firing rate improves the accuracy of the recovery of
the model parameters. Applying DCM to data generated from the original, non-
linear, Jansen-Rit model was less accurate. As with the linear model, reduced
coupling limited the recovery of the parameters of the interneuron
populations. The parameters of a population of uncoupled neurons could be
recovered from a simulated LFP in the absence of spiking dynamics, however,
the introduction of spikes impaired the parameter recovery using DCM. The
error in the recovery of the parameters increased with the frequency of
spikes. The simulated LFP from a network of coupled spiking neurons, with a
network structure reflecting that of the Jansen-Rit model, had a spectrum
similar to that of the transfer function. As parameters in the model were
varied in isolation the variation in the peak amplitude of the spectra was
correlated with that of the transfer function. The frequency at which the peak
occurred was correlated for some, but not all, of the parameters. There was a
similar result for the MLE. Similar to the results for DCM, we found that BMS
could accurately identify which parameter was varied for the linearised
Jansen-Rit model. As follows from the form of the transfer function for the
model, parameters that always appear together as a product could not be
distinguished. The accuracy was less good for data generated from the Jansen-
Rit model without linearisation and from a network of spiking neurons.

Any new methodology requires validation against a standard. Inferring the
properties of the neuronal network that gives rise to an LFP or EEG recording
is challenging. It is difficult to validate results as the properties of the
true network often cannot be confirmed using alternative methods. Modelling
studies have limitations, but allow knowledge of the underling parameters to
evaluate against. Our study investigated the effectiveness of DCM and BMS
applied at steady state, in the frequency domain, using a transfer function
representation. As with previous work, we found DCM and BMS to be effective
when applied to data generated from a linearised Jansen-Rit model. This is to
be expected as, in this form, the model and transfer function are equivalent.
Nonetheless, we note that there are still limitations in this simplified
setting with this particular model as there are groups of parameters that only
appear as a product in the transfer function and so their individual effects
cannot be distinguished. With models that deviate from that relating exactly
to the transfer function, either by incorporating non-linearity into the
response of a Jansen-Rit model or by modelling distinct neurons with spiking
dynamics, we found that DCM and BMS in this form was much less effective. DCM
would be most accurate when applied to networks with a low firing rate and
with as near linear firing rate in response to mean potential as possible. We
found that changes could be identified better than absolute values.

The ability to infer the properties of a network of neurons from a recording
of averaged electrical activity, such as EEG, would be incredibly valuable
both for fundamental neuroscience and clinical investigations. For example,
allowing the identification of a channelopathy causing epilepsy via non-
invasive EEG recording, thereby allowing specific targeting of anti-epileptic
drugs. Our investigation, although limited, found DCM and BMS in the form
studied here to not be effective in achieving this aim. This suggests caution
is required when employing these methods.

Our study was limited in a number of ways. Firstly, as with any modelling-
based study, the synthetic data that we generated had fundamental differences
from true EEG and LFP recordings as our models are extremely simplified in the
dynamics they include. This is, however, not likely to lead to an
underestimate of the effectiveness of the methods as the models used were
specifically formulated to be as similar as possible to the Jansen-Rit model
from which the relevant transfer function was generated. Our approach was
handicapped by employing uninformative flat priors, thereby calculating
maximum likelihood estimates. More informative priors would be likely to
improve accuracy.

Our investigation is specific to the application of DCM and BMS at steady
state, in frequency space, using a transfer function representation of a
Jansen-Rit model. There are other similar approaches and it would be valuable
to investigate their validity using similar modelling studies. It would also
be valuable to test DCM and BMS on models of networks of spiking neurons with
more realistic connectivity and conduction delays between neurons.

Code availability
Code to implement the models presented here can be found in the following
GitHub repository: https://github.com/mgthomas1/Thomas2023

Acknowledgements
I thank Prof. A. Farmery and Prof. R. McShane for their interest and support
in developing this project.

## References

  * Buzsaki et al. 2012 G. Buzsaki, C.A. Anastassiou, and C. Koch.  The origin of extracellular fields and currents - EEG, ECoG, LFP and spikes.  _Nature Reviews Neuroscience_ , 13:407–420, 2012.
  * Friston et al. 2003 K.J. Friston, L. Harrison, and W. Penny.  Dynamic causal modelling.  _NeuroImage_ , 19:1273–1302, 2003.
  * Penny et al. 2004 W.D. Penny, K.E. Stephan, A. Mechelli, and K.J. Friston.  Comparing dynamic causal models.  _NeuroImage_ , 22:1157–1172, 2004.
  * David et al. 2006 O. David, S.J. Kiebel, L.M. Harrison, J. Mattout, J.M. Kilner, and K.J. Friston.  Dynamic causal modelling of evoked responses in EEG and MEG.  _NeuroImage_ , 30:1255–1272, 2006.
  * Moran et al. 2008 R.J. Moran, K.E. Stephan, S.J. Kiebel, N. Rombach, W.T. O’Connor, K.J. Murphy, R.B. Reilly, and K.J. Friston.  Bayesian estimation of synaptic physiology from the spectral responses of neural masses.  _NeuroImage_ , 42:272– 284, 2008.
  * Kiebel et al. 2009 S.J. Kiebel, M.I. Garrido, R. Moran, C-C. Chen, and K.J. Friston.  Dynamic causal modelling for EEG and MEG.  _Human Brain Mapping_ , 30:1866–1876, 2009.
  * Kahan and Foltynie 2013 J. Kahan and T. Foltynie.  Understanding DCM: Ten simple rules for the clinician.  _NeuroImage_ , 83:542–549, 2013.
  * Friston et al. 2019 K.J. Friston, K.H. Preller, C. Mathys, H. Cagnan, J. Heinzle, A. Razi, and P. Zeidman.  Dynamic causal modelling revisited.  _NeuroImage_ , 199:730–740, 2019.
  * Jafarian et al. 2019 A. Jafarian, P. Zeidman, V. Litvak, and K. Friston.  Structure learning in coupled dynamical systems and dynamic causal modelling.  _Philosophical Transactions A_ , 377:20190048, 2019.
  * Moran et al. 2011a R.J. Moran, F. Jung, T. Kumagai, H. Endepols, R. Graf, R.J. Dolan, K.J. Friston, K.E. Stephan, and M. Tittgemeyer.  Dynamic causal models and physiological inference: A validation study using isoflurane anaesthesia in rodents.  _PLoS ONE_ , 6(8):e22790, 2011a.
  * Moran et al. 2009 R.J. Moran, K.E. Stephan, T. Seidenbecher, H.C. Pape, R.J. Dolan, and K.J. Friston.  Dynamic causal models of steady-state responses.  _NeuroImage_ , 44:796–811, 2009.
  * Papadopoulou et al. 2015 M. Papadopoulou, M. Leite, P. van Mierlo, K. Vonck, L. Lemieux, K. Friston, and D. Marinazzo.  Tracking slow modulations in synaptic gain using dynamic causal modelling: Validation in epilepsy.  _NeuroImage_ , 107:117–126, 2015.
  * Cooray et al. 2023 G.K. Cooray, R.E. Rosch, and K.J. Friston.  Modelling cortical network dynamics.  _arXiv_ , arXiv:2301.01144, 2023.
  * Moran et al. 2011b R.J. Moran, M. Symmonds, K.E. Stephan, K.J. Friston, and R.J. Dolan.  An in vivo assay of synaptic function mediating human cognition.  _Current Biology_ , 21:1320–1325, 2011b.
  * Moran et al. 2013 R.J. Moran, P. Campo, M. Symmonds, K.E. Stephan, R.J. Dolan, and K.J. Friston.  Free energy, precision and learning: the role of cholinergic neuromodulation.  _Journal of Neuroscience_ , 33:8227–8236, 2013.
  * Moran et al. 2014 R.J. Moran, M. Symmonds, R.J. Dolan, and K.J. Friston.  The brain ages optimally to model its environment: Evidence from sensory learning over the adult lifespan.  _PLoS Computational Biology_ , 10(1):e1003422, 2014.
  * Symmonds et al. 2018 M. Symmonds, C.H. Moran, M.I. Leite, C. Buckley, S.R. Irani, K.E. Stephan, K.J. Friston, and R.J. Moran.  Ion channels in EEG: Isolating channel dysfunction in NMDA receptor antibody encephalitis.  _Brain_ , 141:1691–1702, 2018.
  * van Wijk et al. 2018 B.C.M. van Wijk, H. Cagnan, V. Litvak, A.A. Kuhn, and K.J. Friston.  Generic dynamic causal modelling: An illustrative application to Parkinsons disease.  _NeuroImage_ , 181:818–830, 2018.
  * Gilbert et al. 2016 J.R. Gilbert, M. Symmonds, M.G. Hanna, R.J. Dolan, K.J. Friston, and R.J. Moran.  Profiling neuronal ion channelopathies with non-invasive brain imaging and dynamic causal models: Case studies of single gene mutations.  _NeuroImage_ , 124:43–53, 2016.
  * Wilson and Cowan 1972 H.R. Wilson and J.D. Cowan.  Excitatory and inhibitory interactions in localized populations of model neurons.  _Biophysical Journal_ , 12:1–24, 1972.
  * Wilson and Cowan 1973 H.R. Wilson and J.D. Cowan.  A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue.  _Kybernetik_ , 13:55–80, 1973.
  * Jansen et al. 1993 B.H. Jansen, G. Zouridakis, and M.E. Brandt.  A neurophysiologically-based mathematical model of flash visual evoked potentials.  _Biological Cybernetics_ , 68:275–284, 1993.
  * Jansen and Rit 1995 B.H. Jansen and V.G. Rit.  Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns.  _Biological Cybernetics_ , 73:357–366, 1995.
  * Ableidinger et al. 2017 M. Ableidinger, E. Buckwar, and H. Hinterleitner.  A stochastic version of the Jansen and Rit neural mass model: Analysis and numerics.  _The Journal of Mathematical Neuroscience_ , 7:8, 2017.
  * Ahmadizadeh et al. 2018 S. Ahmadizadeh, P.J. Karoly, D. Nesic, D.B. Grayden, M.J. Cook, D. Soudry, and D.R. Freestone.  Bifurcation analysis of two coupled Jansen-Rit neural mass models.  _PLoS ONE_ , 13(3):e0192842, 2018.
  * David and Friston 2003 O. David and K.J. Friston.  A neural mass model for MEG/EEG: Coupling and neuronal dynamics.  _NeuroImage_ , 20:1743–1755, 2003.
  * Stephan et al. 2008 K.E. Stephan, L. Kasper, L.M. Harrison, J. Daunizeau, H.E.M. den Ouden, M. Breakspear, and K.J. Friston.  Nonlinear dynamic causal models for fMRI.  _NeuroImage_ , 42:649–66, 2008.
  * FitzHugh 1961 R. FitzHugh.  Impulses and physiological states in theoretical models of nerve membrane.  _Biophysical Journal_ , 1:445–466, 1961.
  * Harrison et al. 2015 P.M. Harrison, L. Badel, M.J. Wall, and M.J.E. Richardson.  Experimentally verified parameter sets for modelling heterogeneous neocortical pyramidal-cell populations.  _PLoS Computational Biology_ , 11(8):e1004165, 2015.

Appendix
Stability analysis

Here we analyse the stability of the linearized Jansen-Rit model. The transfer
function of the Jansen-Rit model considered is:

$H(s)=\left(\frac{(s+a)^{2}}{aA}-\frac{aAC_{1}C_{2}\gamma^{2}}{(s+a)^{2}}+\frac{bBC_{3}C_{4}\gamma^{2}}{(s+b)^{2}}\right)^{-1}$
(16)

If the real parts of all of the poles of the transfer function are negative
the model is stable, if the real part of any of the poles are positive the
system is unstable. While keeping the other parameters constant, we seek the
value of $C$ that results in a transition from all poles having a real part
that is negative to having a pole (or pair of poles) with a real part that is
zero.

We seek $C$ such that there is an $s=0+il$ where the denominator of the
transfer function is zero while the numerator is finite. This yields a pole
with zero real part. For the denominator to be zero both the real and
imaginary components must be zero.

For the real part of the denominator to be zero:

$a^{2}-l^{2}-\frac{a^{2}A^{2}C_{1}C_{2}\gamma^{2}(a^{2}-l^{2})}{(a^{2}-l^{2})^{2}+(2al)^{2}}+\frac{aAbBC_{3}C_{4}\gamma^{2}(b^{2}-l^{2})}{(b^{2}-l^{2})^{2}+(2bl)^{2}}=0$
(17)

For the imaginary part of the denominator to be zero:

$1+\frac{a^{2}A^{2}C_{1}C_{2}\gamma^{2}}{(a^{2}-l^{2})^{2}+(2al)^{2}}+\frac{Ab^{2}BC_{3}C_{4}\gamma^{2}}{(b^{2}-l^{2})^{2}+(2bl)^{2}}=0$
(18)

Noting the relationship: $C=C_{1}=\frac{5}{4}C_{2}=4C_{3}=4C_{4}$, $C$ can be
eliminated by combining equations 17 and 18, yielding a form for $l$:

$(a^{2}-l^{2})(\frac{8}{5}a^{2}A(b^{2}+l^{2})^{2})-\frac{bB}{16}(a^{2}+l^{2})^{2}(b(a^{2}-l^{2})+a(b^{2}-l^{2}))=0$
(19)

This can readily be solved by substituting $y=l^{2}$ and applying the cubic
formula. The solution for $l$ can then be substituted back into the form for
either the real or imaginary part of the denominator to find C. For example:

$C^{2}=\frac{-(a^{2}+l^{2})^{2}(b^{2}+l^{2})^{2}}{0.8(aA\gamma)^{2}(b^{2}+l^{2})^{2}-AB(\frac{b\gamma}{4})^{2}(a^{2}+l^{2})^{2}}$
(20)

Figure 9: Accuracy of the recovery of the parameters of the linearised Jansen-
Rit model using DCM. The relative likelihood of parameter values for
illustrative values of connectivity ($C=50$, black; $C=100$, red; $C=300$,
dashed blue) are shown against the relative parameter value (a parameter
multiplier of one corresponds to the true value). The results here are from
fitting the transfer function to the Fourier transform of the model output
allowing knowledge of the rate of input spikes to the model. Figure 10:
Accuracy of the recovery of the parameters of the linearised Jansen-Rit model
using DCM. The relative likelihood of parameter values for illustrative values
of connectivity ($C=50$, black; $C=100$, red; $C=300$, dashed blue) are shown
against the relative parameter value (a parameter multiplier of one
corresponds to the true value). The results here are from fitting the transfer
function to the Fourier transform of the model output with both the transfer
function and Fourier transform normalised, reflecting the scenario where the
input spike rate is unknown. Figure 11: Accuracy of the recovery of the
parameters of the Jansen-Rit model with the orignal sigmoid function using
DCM. The relative likelihood of parameter values for illustrative values of
connectivity ($C=50$, black; $C=100$, red; $C=300$, dashed blue) are shown
against the relative parameter value (a parameter multiplier of one
corresponds to the true value). The results here are from fitting the transfer
function to the Fourier transform of the model output allowing knowledge of
the rate of input spikes to the model. Figure 12: Accuracy of the recovery of
the parameters of the Jansen-Rit model with the orignal sigmoid function using
DCM. The relative likelihood of parameter values for illustrative values of
connectivity ($C=50$, black; $C=100$, red; $C=300$, dashed blue) are shown
against the relative parameter value (a parameter multiplier of one
corresponds to the true value). The results here are from fitting the transfer
function to the Fourier transform of the model output with both the transfer
function and Fourier transform normalised, reflecting the scenario where the
input spike rate is unknown. Figure 13: Accuracy of the recovery of the
parameters of the linearised Jansen-Rit model (black) and original model (red)
given knowledge of the input spike rate. The relative error in the maximum
likelihood estimates of the parameters is shown as a function of the coupling
constant, $C$, used in the model.  Figure 14: Accuracy of the recovery of the
parameters of the linearised Jansen-Rit model (black) and original model (red)
without knowledge of the input spike rate. To apply DCM both the transfer
function and Fourier transform were normalised by their means. The relative
error in the maximum likelihood estimates of the parameters is shown as a
function of the coupling constant, $C$, used in the model. Figure 15: Recovery
of the parameters of networks of spiking neurons. Each row corresponds to
variation of a single parameter as indicated ($a$, $A$, $b$, $B$, $C_{3}$, and
$C_{4}$). The first column shows the transfer functions for the equivalent
linearised Jansen-Rit model as the relevant parameter is varied from 50% (red)
to 200% (magenta). The second column shows the Fourier transform of the model
LFP as the relevant parameter is varied from 50% (red) to 200% (magenta). The
first and second columns have a logarithmic scale on the y-axis. Subsequent
columns show the relationship between the frequency at which the peak in the
spectra and the peak in the transfer function occurs, the relative amplitude
of the peak in the spectra and the transfer function, and the true parameter
and the maximum likelihood estimate. In each case the solid red line shows a
least squares fit to the data points and the dashed black line shows the line
$y=x$.