# TVB C++: A Fast and Flexible Back-End for The Virtual Brain

Ignacio Martín  ViRVIG, Universitat de Girona, Girona, Spain Gorka Zamora-
López  Center for Brain and Cognition, Computational Neuroscience Group,
Department of Information and Communication Technologies, Universitat Pompeu
Fabra, Spain Jan Fousek  Institut de Neurosciences des Systèmes, Aix Marseille
Université, Marseille, France Michael Schirner  Berlin Institute of Health,
Charité – Universitätsmedizin Berlin, Germany Petra Ritter  Berlin Institute
of Health, Charité – Universitätsmedizin Berlin, Germany Viktor Jirsa
Institut de Neurosciences des Systèmes, Aix Marseille Université, Marseille,
France Gustavo Deco  Center for Brain and Cognition, Computational
Neuroscience Group, Department of Information and Communication Technologies,
Universitat Pompeu Fabra, Spain Institució Catalana de la Recerca i Estudis
Avançats (ICREA), Barcelona, Spain Gustavo Patow  ViRVIG, Universitat de
Girona, Girona, Spain Center for Brain and Cognition, Computational
Neuroscience Group, Department of Information and Communication Technologies,
Universitat Pompeu Fabra, Spain e-mail<EMAIL_ADDRESS>
###### Abstract

This paper introduces TVB C++, a streamlined and _fast_ _C++_ Back-End for The
Virtual Brain (TVB), a renowned platform and a benchmark tool for full-brain
simulation. TVB C++ is engineered with speed as a primary focus while
retaining the flexibility and ease of use characteristic of the original TVB
platform. Positioned as a complementary tool, TVB serves as a prototyping
platform, whereas TVB C++ becomes indispensable when performance is paramount,
particularly for large-scale simulations and leveraging advanced computation
facilities like supercomputers. Developed as a TVB-compatible Back-End, TVB
C++ seamlessly integrates with the original TVB implementation, facilitating
effortless usage. Users can easily configure TVB C++ to execute the _same_
code as in TVB but with enhanced performance and parallelism capabilities.

###### keywords:

The Virtual Brain, Whole-Brain model, C++, Fast

##  1 Introduction

In recent years, technological advancements in noninvasive neuroimaging have
opened unprecedented windows into the inner workings of the human brain.
Techniques such as EEG, MEG, fMRI, and PET have revolutionized our
understanding of perception, cognition, and behavior, enabling detailed
analyses from the microscopic to the macroscopic scale.

Within this landscape, computational tools have played a pivotal role in
unraveling the complexities of neural dynamics. While various frameworks
exist, including those focusing on microscopic neuron models [1, 2, 3, 4, 5]
or mesoscopic systems [6, 7], the emphasis on whole-brain modeling has led to
the emergence of platforms like The Virtual Brain (TVB) [8, 9, 10, 11, 12,
13]. TVB has garnered widespread recognition for its ability to simulate,
analyze, and infer neurophysiological mechanisms across different brain
scales.

Despite its great versatility, TVB faces a competitive challenge in terms of
computational speed when compared to other simulation platforms, especially
when dealing with extensive networks or large-scale simulations. Several
approaches have been proposed in recent years to boost the performance of TVB,
e.g., by porting specific models into the C-language [CITE] or interfacing
with CUDA [CITE]. However, these attempts lacked the flexibility of TVB
because they were restricted to specific scientific questions they wanted to
address, without resulting in a general setting with the flexibility TVB has.

Here, we introduce TVB C++, a groundbreaking addition to the TVB family aimed
at addressing these limitations simultaneously. Developed upon the powerful
C++ programming language and the advanced _Eigen_ numeric library [14], TVB
C++ offers great computational performance without compromising on flexibility
or ease of use. Unlike previous efforts —focused solely on performance—TVB C++
takes a holistic approach by retaining the same structure of the original TVB
platform. This allows for TVB C++ to seamlessly integrate with TVB through
Python bindings and, therefore, effortlessly complement existing workflows,
offering a swift and efficient backend solution for bulk simulations where
computational time is critical. With TVB C++, researchers gain access to a
cutting-edge tool that empowers them to explore the intricacies of large-scale
brain dynamics or conduct extensive parameter explorations, TVB C++ stands
ready to help advance the field, offering unprecedented speed, performance,
and versatility.

The paper is organized as follows. Section 2 provides a brief overview of
existing computational tools to simulate neural and brain activity at
different scales. In Section 3 we present the original TVB implementation, its
purpose, and general structure, and in Section 4 we present TVB C++ in detail,
with all its input data formats, code structures, models, observables, and,
finally, its Python bindings. In Section 5 we present some use cases for TVB
C++, along with the respective results obtained, and finally we summarize the
impact and possibilities of TVB C++ in Section 6.

##  2 Background

In recent decades, advances in neuroimaging have facilitated unprecedented
levels of analysis across various scales, ranging from the microscopic realm,
which entails simulating the complex interplay among ionic channels,
neurotransmitters, receptors, and action potentials, to the macroscopic scale
encompassing whole-brain simulations.

Comprehensive investigations enable the mechanistic comprehension of
population dynamics, contributing to the elucidation of complex structures
such as Resting State Networks [15, 16, 17], cognition [18, 19], and
pathologies like Epilepsy [20, 21], Stroke [22, 23, 24], or Alzheimer’s
Disease [25, 26]. It is within the realm of whole-brain simulations that true
integration occurs, leading to the emergence of cognition. Consequently, these
models offer a panoramic view of neuronal dynamics at the global macroscopic
level. Computational tools have evolved to simulate all these effects, both at
the algorithmic and at the computational power levels. This has enabled the
simulation of large neuronal models that, in turn, enabled accurate
predictions illuminating both theoretical and clinical aspects of our brains.
Among the different software frameworks developed to this end, we can
distinguish three big groups:

  * •

Class I, or microscopic-level simulators, which aim at simulating the activity
of individual neurons with high detail, such as NEST [1], Brian [2], NetPyNE
[3], and NEURON [4, 5]. These frameworks have proven particularly useful for
simulating large networks of spiking neurons. In particular, NEST and Brian
had their main target on point neurons, while NetPyNE and NEURON focused
primarily on morphologically extended neurons. These frameworks are usually
used to simulate and predict the behavior of large networks of neurons, from
which neural mass models are often derived.

  * •

Class II, or mesoscopic-level simulators, intends to simulate networks of
interconnected neurons to simulate microscopic and mesoscopic neural circuits
such as the internals of subcortical ganglia or cortical columns. There are
other frameworks created with the specific objective of simulating mesoscopic
systems, such as NENGO [6], focusing on applications in cognitive science, or
Brain Modeling Toolkit [7], focusing on multiscale neural population circuits.
Despite their obvious and many advantages, these software frameworks are most
specifically designed for whole-brain modeling, so they are rarely used in
this context given the difficult process of calibrating the models.

  * •

Class III, or whole-brain simulators, aims at simulating whole-brain network
activity at the resolution of interconnected brain regions. The main
ingredients of these models are that a neural population model simulates each
brain region (ROI), typically interconnected via the white matter fibers as
captured by tract tracing or tractography based on diffusion imaging. One of
the first and main open platforms for this endeavor is _The Virtual Brain_
(TVB)111https://www.thevirtualbrain.org/, which, since its inception, has
gained wide recognition as the main platform for full-brain simulations, being
a key element in the world’s first clinical trial on predictive brain modeling
in epilepsy surgery [9, 8]. Other options are PyRates [27], and neurolib [28],
which have been successfully applied to simulate a few use cases.

The Virtual Brain is the most extended and popular framework for Class III
whole-brain modeling. It is, in essence, a neuroinformatics platform with a
brain simulator that incorporates a range of neuronal models and dynamics at
its core. It seamlessly integrates computational modeling and multimodal
neuroimaging tools [9, 8], allowing the model-based simulation, analysis, and
inference of neurophysiological mechanisms over several brain scales that
underlie the generation of macroscopic neuroimaging signals. With TVB,
creating personalized virtual brains is possible [29, 10], as well as studying
the intricate multi-scale neural mechanisms within the brain [30]. The most
recent advances in TVB include its integration with cloud services through the
European platform EBRAINS [11].

However, the original TVB implementation —which we will refer to as _TVB-
Python_ in the following— has in computational speed its biggest drawback. To
overcome slow execution speeds, two different strategies have been
implemented. The first one, TVB-
HPC222https://wiki.ebrains.eu/bin/view/Collabs/rateml-tvb [12], automatically
produces high-performance codes for CPUs and GPUs using an easy XML-based
language called RateML for model specification. RateML is based on the domain-
independent language _LEMS_ [31], which allows for the declarative description
of computational models using a simple XML syntax. The already existing
example implementations can be easily adapted to test different models,
without requiring any knowledge about algorithmic optimization. The second
one, Fast_TVB (first used in [30]), is a specialized high-performance
implementation of the "Reduced Wong Wang" model [32] written in C. It makes
use of several optimization strategies and a sparse memory layout to
efficiently use CPU resources, which makes it possible to simulate extremely
large models with millions of nodes even on a standard computer in a
reasonable time. Aquilué-Llorens et al. [33] presented a CPU-based version of
the TVB-AdEx model [34, 35], used to study mechanisms underlying the emergence
of conscious- and unconscious-like brain state dynamics using an MPI multi-
node CPU setup. More recently, [13] presented the TVB-HPC framework, a modular
set of methods used here to specifically implement the TVB-AdEx model for GPU
and analyze emergent dynamics, accelerating simulations and thus reducing
computational resource requirements, which in turn enabled larger parameter-
space explorations.

However, despite its broad acceptance as _the_ reference platform for whole-
brain simulations, it still exhibits some limitations when it comes to very
large brain simulations, mainly in terms of its speed. For instance, a
simulation with a parcellation with around 400 ROIs using the Balanced
Excitation-Inhibition (BEI) model [32] using BOLD signals of around 10-minute
recordings (approximately XX temporal samples, or volumes) takes around 5
minutes in a high-end processor. Suppose this evaluation is inserted into a
fitting procedure where hundreds of evaluations take place. In that case, the
total time for a single experiment on a single subject may rise to 1 week,
which can become a serious problem even if using powerful parallel
supercomputers.

##  3 TVB Python

The Virtual Brain (TVB) stands as a neuroinformatics software platform
designed for whole-brain network simulations, leveraging biologically
realistic connectivity to facilitate model-based inference of
neurophysiological mechanisms across diverse brain scales. See Figure 1.
Through the incorporation of such connectivity, TVB enables the generation of
macroscopic neuroimaging signals, including functional MRI (fMRI), EEG, and
MEG, as depicted in Figure 1. By harnessing individual subject data, TVB
empowers the reproduction and assessment of personalized brain configurations.
This personalized approach fosters exploration into the ramifications of
pathological alterations within the system, thereby facilitating
investigations into potential strategies to mitigate adverse processes. Within
TVB, a biologically plausible, large-scale connectivity of brain regions is
integrated. This connectivity is mediated through long-range neural fiber
tracts identified via tractography-based methods. One differential trait of
TVB and TVB C++, is their capability to handle connectivity both at the
strength level, and also the time-delays associated to fiber lengths. In
practice, the implementation for weights and tract lengths entails full
matrices without any inherent symmetry constraints.

From a high-level perspective, the architecture of TVB, as illustrated in
Figure 2, accommodates interaction through various interfaces tailored to
specific user needs and deployment scenarios. Structurally, TVB comprises two
primary components: the backend, referred to as the "scientific computing
core," and the supporting framework, which includes the graphical user
interface and data visualization tools. TVB relies on an efficient and
resilient storage data system to manage both metadata information and the
underlying data. Functionally, TVB backends are designed to seamlessly support
various interfaces, ensuring consistent functionality regardless of the
selected interface.

In essence, a TVB application defines the logic necessary for data input,
processing based on configurable parameters, and visualization of results. In
the original Python implementation of TVB, while a single simulation step may
not be readily distributable across multiple processing threads, each
simulation step can be computed independently in different threads,
potentially across distinct nodes. TVB documentation advises the deployment of
complex simulation tasks on clusters to mitigate computational expenses.

As previously mentioned, TVB-Datatypes serve as a middleware layer
facilitating the management and flow of data between the scientific kernel and
the supporting framework. Essentially, TVB-Datatypes are annotated data
structures encompassing one or more data attributes alongside associated
descriptive information, inclusive of requisite methods for data manipulation.

The simulation core, the final module in our delineation, utilizes structural
connectivity information within a comprehensive model of neural dynamics
across the entire brain. Here, users can delineate the spatial structure upon
which the model operates, as well as the hierarchical connections between
nodes through the structural connectivity framework. Subsequently, the system
numerically integrates the resultant coupled system of differential equations
to replicate emergent brain dynamics. Central to this process, the Simulator
class serves as both a repository for all pertinent information and as the
conductor orchestrating the various steps associated with the simulation, all
transparently to the user.

##  4 TVB C++

In response to the constrained performance of TVB Python, particularly evident
when handling extensive networks or simulations necessitating extensive
parameter exploration, TVB C++ emerges as a simulator tailored for whole-brain
network models (BNMs).

Reproducing the functionality of TVB-Python, TVB C++ can simulate neural
population models coupled by structural connectivity, which is an aggregated
representation of the brain’s white matter axon fiber bundle network. In
addition to being implemented as a C++ library, TVB C++ offers Python bindings
to seamlessly plug it into existing pipelines, and subsequently serve as a
fully integrated TVB Back-End. This way, by simply configuring the TVB
simulator to use TVB C++ instead of the default implementation in Python,
existing code can be used without changes. The general architecture of TVB C++
as a back-end for TVB is illustrated in Figure 3.

The architecture of TVB C++ centers on two distinct interfaces tailored to
user interaction. Currently, the console interface operates independently of
the storage layer, storing results solely in memory. Consequently, users
utilizing the console interface must manually oversee data import and export
operations. On the other hand, its back-end comprises blocks that seamlessly
interface with various top application layers. The _data types_ serve as the
common language between different components such as analyzers, visualizers,
simulators, and uploaders. These datatypes embody "active data," signifying
that when TVB C++ is configured with a database, instances encapsulating data
types are automatically retained.

In the following subsections, we will introduce the main TVB C++ components,
namely its inputs, code structure, whole-brain models, integrators, history,
coupling, and monitors. All these components represent functional modules that
are completely user-configurable (see Section 3), allowing the users to
interact with them to tailor their simulation and analysis needs.

###  4.1 Inputs

In general, TVB C++ can load brain model data from a variety of sources and
formats, with easy expandability to new formats and registration modalities.
One of the main inputs for brain data consists of a structural connectome SC
(in our case, weights files) and a parameter set (parameter file), plus
spontaneous blood-oxygenation-level-dependent (BOLD) as stored in functional
MRI data (fMRI), together with some regional information that depends on the
problem at hand, from neuroreceptor data to gene expressions. Optionally, the
user can provide information about the lengths of the tracts between the
nodes, which would introduce _delays_ in the communications between the
regions. This is also considered in the TVB C++ implementation.

####  4.1.1 Data formats

While TVB allows loading data from Matlab data files (MAT), our version uses
the Numpy (NPZ) format. We provide a Python script to convert MAT files into
NPZ files.

###  4.2 Code structure

TVB C++ replicates the same structure of the original TVB library. That means
that the structure of the classes is the same, providing the same
functionality and flexibility. However, TVB C++ only focuses on the simulation
part of the original TVB, and it does not provide any user interface. The two
main parts of the library are the simulator module modulename.py, consisting
of the population models to simulate the activity of one brain region, e.g.,
mean-field AdEx, Wong-Wang, FRE; coupling, which can be either dense or
sparse; the integrator, where the different integration schemes are
implemented; the monitors, including raw data monitor and the temporal
averaging monitor, as well as a convolution-based BOLD; the noise
implementation; the BOLD conversion system based on the Ballon-Windkessel
model; and the simulator module itself, in charge of wrapping all the other
components in a single, uniform shell. On the other hand, the other big
module, the tools, is composed of a series of submodules, being the
observables, e.g., FC, FCD, and the Phase Interaction Matrices a few examples,
as well as the interfaces with numpy and Scipy. See Figure 4.

While TVB is heavily based on NumPy, no version of this library exists in C++.
To provide similar functionality we have chosen the _Eigen_ library [14],
specifically version 3.9 since it provides advanced matrix indexing similar to
NumPy. Moreover, Eigen takes advantage of SIMD instruction sets, and it can be
combined with BLAS/LAPACK, Intel MKL library, and CUDA.

All the objects that can be configured in a simulation are:

  * •

Whole-Brain Models

  * •

Integrators

  * •

History

  * •

Coupling

  * •

Monitors

###  4.3 Population and mean-field models of regional activity

TVB C++ has a modular design that allows the incorporation of a large variety
of whole-brain models. Model definitions are very similar to TVB, and with TVB
C++ this implies deriving from the abstract class tvb::Model. The main method
to implement is tvb::Model::operator() (equivalent to the dfun() method in
TVB) that computes the derivatives of the state variables.

Currently, several methods have been implemented and tested, namely, the
Balanced Excitation-Inhibition (BEI) model by Deco et al. [32], the Montbrió,
Pazó and Roxin’s FRE model [36, 37], and the mean-field AdEx [34, 38, 39, 40].
See Section 5.

###  4.4 Numerical integrators

TVB allows the use of different integrators to compute simulations with the
same model. TVB C++ follows the same approach and provides an abstract class
Integrator that can be derived and this way implements different integration
strategies. Right now we provide the Euler deterministic and stochastic
(Euler-Maruyama) integrators. In this latter case, stochastic integrators can
be combined with different noise generators, in the same way, that TVB does.

###  4.5 History and time-delays

TVB takes into account the delay in signal transmission between nodes of the
brain model. The delays are computed from the tract length matrix and the
signal speed. TVB C++ also keeps track of the signal history and provides two
modes: a dense history that keeps track of all signals at all possible delays,
and a no-delay history, which is more efficient and can be used in some
models.

###  4.6 Coupling

The coupling functions allow activity (state variables) that have been
propagated over the long-range connectivity to pass through these functions
before entering the equations (tvb::Model::operator()) describing the local
dynamics.

At the moment, only linear coupling is implemented, but other functions are
straightforward to use.

###  4.7 Monitors

Another feature of TVB is the Monitors, objects that record significant values
from the simulation. In their simplest form, they return all the simulated
data, subsampled data, spatially averaged temporally subsampled, or temporally
averaged subsamples. The current version of TVB C++ implements the two former
ones.

TVB also has more elaborate monitors that instantiate a physically realistic
measurement process on the simulation, such as EEG, MEG, and fMRI (BOLD). In
TVB C++ we follow a slightly different approach and these measurement
processes are performed after the simulation, using the monitor data as input.

###  4.8 Extra tools

To make TVB C++ a useful tool, the ability to simulate Whole-Brain models is
just the first step. There are many extra tools needed to analyze the results
of such simulations. TVB C++ provides some of these tools, and provides an
example implementing all the necessary steps to reproduce the main results of
the corresponding models enumerated above and described in Section 5: The BEI
model [32], the FRE model [36, 37] and the mean-field AdEx model [34, 38, 39,
40]. We have chosen this work since it allows us to compare standard features
already existing in TVB and high-performance implementations. To reproduce
these works, some steps need to be taken:

  * •

Perform simulations using the corresponding model

  * •

Estimate the BOLD signal from the simulated physiological signals.

  * •

Compute some sort of observable, e.g., sliding-window correlations, from the
obtained BOLD signals. See Section 4.8.2.

  * •

Use a measure, e.g., the $L^{2}$ metric, between the empirical and simulated
observables obtained.

For the last two items, we had to implement some of the features of NumPy and
SciPy into TVB C++, like the Kolmogorov statistic.

####  4.8.1 Estimate BOLD signals from neural output

Once we have obtained the simulated mean field activity, we need to transform
it into a Blood Oxygenation Level Dependent (BOLD) signal. Given that the
conversion from electrical activity to a BOLD signal is, in its mathematical
sense, a low-pass filter, this process can be computed in two different ways:
as a convolution with Volterra kernels, or as a generalized hemodynamic model
such as the Balloon-Windkessel model. In the following, we describe both
implementations, as TVB C++ implements the two methods to create a BOLD signal
from a raw signal produced by the simulation, as the original TVB
implementation also does.

##### Haemodynamic Response Function Model

One option is to use convolutions with Hemodynamic Response Function (HRF) to
compute the output BOLD signal, derived using generic techniques from
nonlinear system identification [41, 42, 43, 44, 45]. The approach employed
allows for the estimation of Volterra kernels, which elucidate the
relationship between stimulus presentation and consequent hemodynamic
responses. This technique represents the BOLD signal as the outcome of a
convolution with a Volterra series [46], which are, essentially, Taylor
expansions extended to encompass dynamical input-state-output systems and
consider the impact of the input at all times in the recent past [42]. As a
result, they are regarded as model-independent, capable of representing any
such dynamical system. Serving as high-order extensions of linear convolution
or "smoothing," Volterra kernels provide a nonlinear characterization of the
hemodynamic response function. Consequently, they facilitate the modeling of
responses to stimuli across various contexts (e.g., different rates of word
presentation) and interactions among stimuli. In TVB C++ this is implemented
as a special case of a _Monitor_ class, which are tools to record significant
values from the simulation, in this case using HRF kernels. As we are
following TVB’s main implementation, we are also using the first-order
Volterra kernel, too.

##### Balloon-Windkessel Model:

We implemented the generalized hemodynamic model of Stephan et al. [47, 48].
We compute the BOLD signal in the $k$-th brain area from the firing rate of
the excitatory pools $H^{(E)}$, such that an increase in the firing rate
causes an increase in a vasodilatory signal, $s_{k}$, that is subject to auto-
regulatory feedback. Blood inflow $f_{k}$ responds in proportion to this
signal inducing changes in blood volume $v_{k}$ and deoxyhemoglobin content
$q_{k}$. The equations relating to these biophysical variables are:

$\begin{split}\frac{ds_{k}}{dt}&=0.5r_{k}^{(E)}+3-ks_{k}-\gamma(f_{k}-1)\\\
\frac{df_{k}}{dt}&=s_{k}\\\
\tau\frac{dv_{k}}{dt}&=f_{k}-v_{k}^{\alpha^{-1}}\\\
\tau\frac{dq_{k}}{dt}&=f_{k}\frac{1-(1-\rho)^{f_{k}^{-1}}}{\rho}-q_{k}\frac{v_{k}^{\alpha^{-1}}}{v_{k}}\end{split}$
(1)

with finally

$B_{k}=v_{0}\left[k_{1}(1-q_{k})+k_{2}(1-\frac{q_{k}}{v_{k}})+k_{3}(1-v_{k})\right]$

being the final measured BOLD signal.

We used the updated version [49] described later, which consists of
introducing the change of variables $\hat{z}=lnz$, which induces the following
change for $z=f_{k}$, $v_{k}$ and $q_{k}$, with its corresponding state
equation $\frac{dz}{dt}=F(z)$, as:

$\frac{d\hat{z}}{dt}=\frac{d\;ln(z)}{dz}\frac{dz}{dt}=\frac{F(z)}{z}$

which results in $z(t)=exp(\hat{z}(t))$ always being positive, which
guarantees proper support for these non-negative states, and thus numerical
stability when evaluating the state equations during evaluation.

####  4.8.2 Observables

##### Functional Connectivity:

It is computed as the Pearson correlation between the model and empirical
estimates of static (i.e., time-averaged) Functional Connectivity, estimated
across all pairs of brain regions (FC); more specifically, the correlation
between the values in the upper triangles of the model and empirical FC
matrices.

##### Sliding-Window Functional Connectivity Dynamics (swFCD):

In [50], to take into account the spatiotemporal fluctuations in functional
brain dynamics over time, the model is fitted to the spatiotemporal dynamics
of the data (i.e., to the Functional Connectivity Dynamics [FCD]) [17, 16,
51]. The measure used is a sliding-window correlation between BOLD signals. To
implement this, TVB C++ re-implements some functions from NumPy and SciPy.

##### Phase Functional Connectivity Dynamics (phFCD):

From each BOLD time series of each ROI, we can Hilbert-transform them to yield
the phase evolution of the regional signals, which can be averaged, to finally
be used to compute the KS distance from other identically processed series.
This yields the Phase Functional Connectivity Dynamics (phFCD), which is a 3D
matrix of $N\times N\times T$-size, where $N$ is the number of nodes in the
parcellation used and $T$ indicates the number of image volumes acquired in
all considered sessions.

###  4.9 Python bindings

Compiling and executing C++ code can be challenging to most neuroscientists
who are more used to user-friendly interfaces. To provide an easy entry point
to use the power of TVB C++, we have developed a Python package that provides
easy access to the features of this library. This package provides bindings
for the different library components, allowing the selection of different
models (e.g., FRE, Dynamic Mean Field, mean-field AdEx), different monitors
(e.g., two different BOLD monitors), integrators (e.g., Euler or Euler-
Maruyama) and observables (e.g., FC, swFCD, phFCD). The API supports
interaction with the original TVB simulator python class so that it can be
used as an alternative Back-End without the need for code adaptation.

The library is hosted in a Github public repository [52] and there is a folder
that contains the Python binding including a _Jupyter notebook_
333https://github.com/neich/tvb-root-cpp/blob/main/python/tvbcpp.ipynb that
makes it easier to test the library. The API as of now is straightforward and
it allows to select of different neural population models, monitors,
integrators, and observables, following the protocols and standards of the
original TVB implementation, thus allowing for seamless integration as a
backend for TVB-python.

In Listing 1, a comprehensive example of TVB C++ Python bindings is provided.
It is essential to note that the TVB C++ module operates as a state machine,
requiring the configuration of all state variables before execution, starting
with the integration variables and loading the structural connectivity matrix
(lines 6-14). Subsequently, a neuronal mass model, FRE model with 6 state
variables, is configured, with noise added only to the last two via the Euler-
Maruyama (referred to as _EulerStochastic_) integrator (lines 16-21). The
initial state of the model’s variables is set to zero (lines 22-25), followed
by the configuration of global coupling (line 27). While the current
implementation only offers linear coupling, the library is designed for easy
extension. Users can query model parameters before execution by utilizing the
get_model_parameters() command (line 29), and change them before execution. An
example of this can be seen below, in Section 5.3, where an example of a
change of some complex model parameters is presented.

A monitor is then added to process the raw simulation data, averaging over 1
ms of state variable 0 (line 32). Finally, the simulation is run with the
run_sim command, specifying the time interval in milliseconds (line 34). The
simulation output comprises a list of monitor results, with each result being
a tuple containing two arrays: a 1D array with time samples and a 2D array
with simulation output data, such as time-series data.

To convert the simulation output to a BOLD signal, three methods are
available: Stephan2007, Stephan2007b, and the original TVB convolution-based
monitor (lines 41-43). Users can select their desired monitor by changing the
appropriate keyword.

Finally, the resulting signal is subsampled for illustrative purposes, and the
results can be plotted using the code snippet provided in Listing 2.

##  5 Use cases

In this section we provide practical use cases implemented within TVB C++.
First, we present a modified version of the Balanced Excitation-Inhibition
model [32], followed by the implementation of the Mean Field Adaptive
Integrate and Fire Neurons (mean-field AdEx) in its first order [35] and
second order [34] variants. Finally, the FRE model [36, 37] is described.

We have tested the TVB C++ version against TVB-Python version 2.8.1, the
FAST_TVB implementation, and RateML. The comparison involves two main tasks,
the ones involving whole-brain simulations using a neural model, and other
computations that use mainly external functions, like BOLD signal generation.
It must be taken into account that FAST_TVB only offers a reduced set of the
full TVB functionality, so we only compared the functionality that was
possible for us to compare. All comparisons were made in a computed with an
Intel i9-13900H CPU (14 cores) and an NVidia laptop RTX 4070 GPU.

###  5.1 Modified Balanced Excitation-Inhibition

In this section, we will use TVB C++ to reproduce the results previously
obtained by Deco et al. [50]. The full Balanced Excitation/Inhibition (BEI)
model [32] at the whole-brain level is expressed by the following system of
coupled differential equations:

$\begin{split}I_{k}^{(E)}&=W_{E}\,I_{o}+w_{+}\,J_{N}\,S_{k}^{(E)}+J_{N}G\sum_{j}C_{kj}S_{j}^{(E)}-J_{k}S_{k}^{(I)}+I_{ext},\\\
I_{k}^{(I)}&=W_{I}\,I_{o}+J_{N}S_{k}^{(E)}-S_{k}^{(I)}+\lambda
J_{N}G\sum_{j}C_{kj}S_{j}^{(E)},\\\
r_{k}^{(E)}&=H^{(E)}(I_{k}^{(E)})=\dfrac{M_{k}^{E}(a_{E}I_{k}^{(E)}-b_{E})}{1-\exp(-d_{E}M_{k}^{E}(a_{E}I_{k}^{(E)}-b_{E}))},\\\
r_{k}^{(I)}&=H^{(I)}(I_{k}^{(I)})=\dfrac{M_{k}^{I}(a_{I}I_{k}^{(I)}-b_{I})}{1-\exp(-d_{I}M_{k}^{I}(a_{I}I_{k}^{(I)}-b_{I}))},\\\
\dot{S}_{k}^{(E)}&=-\dfrac{S_{k}^{(E)}}{\tau_{E}}+(1-S_{k}^{(E)})\,\gamma
H^{(E)}(I_{k}^{(E)}),\\\
\dot{S}_{k}^{(I)}&=-\dfrac{S_{k}^{(I)}}{\tau_{I}}+H^{(I)}(I_{k}^{(I)})\end{split}$

Here, the last _two_ equations should add, when integrating, an uncorrelated
standard Gaussian noise term with an amplitude of $\sigma=0.01nA$ (Euler-
Maruyama integration).

The parameters of this model were originally defined in Wong and Wang [53] to
emulate resting-state conditions, such that each isolated node displays the
typical noisy spontaneous activity with low firing rate ($H^{(E)}\sim 3Hz$)
observed in electrophysiology experiments. Moreover, following Deco et al.
[32], the inhibition weight, $J_{n}$, was adjusted separately for each node
$n$ such that the firing rate of the excitatory pools $H^{(E)}$ remains
clamped at 3Hz even when receiving excitatory input from connected areas. It
has been demonstrated that this mechanism leads to a better prediction of the
resting-state FC and a more realistic evoked activity [32]. We refer to this
model as the balanced excitation-inhibition (BEI) model.

To study the effects of LSD, Deco and co-authors [50] used local regional
information, introduced in the form of a detailed map of $5-HT_{2A}$ receptor
density of the neuromodulator serotonin, obtained from a new high-resolution
human brain in vivo atlas. For this, they defined a global gain scaling
parameter, added to the original fixed gain parameters (which were the same
for all regions), and thus scaling the regional $5-HT_{2A}$ receptor values
influencing the recursive circuits of excitatory and inhibitory neurons. In
mathematical terms, we modify the value of $M_{k}^{E}$ in the above equations
as

$M_{k}^{E}=g_{E}(1+s\;\;5HT2A_{k})$

with $g_{E}$ the value previously used [32], $5HT2A_{k}$ the neuroreceptor
density for the $k$-th region, and $s$ the scaling factor. The procedure
follows the somewhat standard pipeline in these cases: first, fit the model to
the placebo condition but not the LSD condition, i.e., assuming zero values of
the scaling that correspond to the original gain values. Then, a scaling value
that fits the LSD condition (using the sensitivity of the functional
connectivity dynamics) while still using the original whole-brain placebo
model, but now including the new element of receptor binding, which modulates
each region with the different empirical measures of the receptor binding.
Thus, here we verify their results with TVB C++, showing that neurotransmitter
modulation of whole-brain activity dynamics can be quantitatively ascribed to
one type of receptor binding (here $5-HT_{2A}$) that modulates brain-wide
neural responses.

As a first test, we compared TVB C++ with TVB-Python version 2.8.1, with
FAST_TVB, and with RateML. In all cases, we use a reduced Wong Wang model [32]
and the Julich-Brain probabilistic cytoarchitectonic maps as a microstructural
reference parcellation from the multilevel Human Brain Atlas, of 379 regions.
All simulations evaluate 10 seconds and use a deterministic integration with
an integration delta of 0.1 ms. Execution times are shown in Table 1. In
Figure 5, top row, we show the plot of the first state variables ($S_{e}$) of
the 379-region model over 10 seconds. Both results shown are numerically
identical. For the sake of completeness, although the random number generators
in Python and C++ are different, we compared one of the models, the BEI model,
with noise, both for TVB and TVB C++, leading to practically indistinguishable
results except for the noise itself, see Figure 6.

TVB | FAST_TVB | TVB C++
---|---|---
40.1 | 5.1 | 6.4
Table 1: Execution times for the reduced Wong Wang model simulation (in
seconds).

##### Feedback Inhibition Control:

As a use case beyond brain simulations, we have implemented the Feedback
Inhibition Control mechanism proposed by Deco et al. [32]. We have implemented
the optimization method that allows finding the $J_{i}$ values that allow the
neuron’s firing rates to stay in a given interval, for every value of the
global coupling value $G$. In particular, we have reproduced Figure 2C of Deco
et al. [32], see Figure 5, lower right. In this case, the results are not
numerically identical since we are not using the same data as used in the
paper for the test, but, as can be seen, the outcome is completely equivalent.

###  5.2 Mean Field Adaptive Integrate and Fire Neurons (mean-field AdEx)

In this section, we will elaborate on the first and second-order Mean-Field
Adaptive Integrate and Fire (mean-field AdEx) model.

##### First Order model:

A mean-field model of conductance-based networks of adaptive exponential
integrate-and-fire neurons [35]

$\displaystyle T\frac{\partial\nu_{\mu}}{\partial t}=$ $\displaystyle-
F_{\mu}(\nu_{e},\nu_{i})+\nu_{\mu},{\forall}\mu\in\\{e,i\\}$
$\displaystyle\frac{\partial W_{k}}{\partial t}=$ $\displaystyle
W_{k}/tau_{w}-b*E_{k}$

where $\nu_{\mu}=\\{E,I\\}$, and $F_{\lambda}$ is the Transfer Function, which
should be tailored for each specific type of single neuron model desired.
Carlu et al. [39] used this framework to obtain Transfer Functions for diverse
models such as nonlinear integrate-and-fire or Hodgkin-Huxley models.

##### Second Order model:

A mean-field model of conductance-based networks of spiking neurons with
adaptation [34]. The general formulation for the mean-field AdEx adaptation
second order model (a.k.a., Zerlaut’s model, di Volo’s model, …) as a
dynamical unit at a node $k$ in a network with $l$ nodes reads
$\forall\mu,\lambda,\eta\in\\{e,i\\}^{3}\,$:

$\displaystyle T\,\frac{\partial\nu_{\mu}}{\partial t}=$
$\displaystyle(\mathcal{F}_{\mu}-\nu_{\mu})+\frac{1}{2}\,c_{\lambda\eta}\,\frac{\partial^{2}\mathcal{F}_{\mu}}{\partial\nu_{\lambda}\partial\nu_{\eta}}$
$\displaystyle T\,\frac{\partial c_{\lambda\eta}}{\partial t}=$ $\displaystyle
A_{\lambda\eta}+(\mathcal{F}_{\lambda}-\nu_{\lambda})\,(\mathcal{F}_{\eta}-\nu_{\eta})+$
$\displaystyle
c_{\lambda\mu}\frac{\partial\mathcal{F}_{\mu}}{\partial\nu_{\lambda}}+c_{\mu\eta}\frac{\partial\mathcal{F}_{\mu}}{\partial\nu_{\eta}}-2c_{\lambda\eta}$
$\displaystyle\frac{\partial W_{k}}{\partial t}=$ $\displaystyle
W_{k}/tau_{w}-b*E_{k}$

with:

$A_{\lambda\eta}=\begin{cases}\frac{\mathcal{F}_{\lambda}\,(1/T-\mathcal{F}_{\lambda})}{N_{\lambda}}\qquad&\textrm{if
}\lambda=\eta\\\ 0\qquad&\textrm{otherwise}\end{cases}$

This model has recently been extended to Whole-Brain modeling [40] by adding
the coupling through the structural connectivity into the synaptic current.

We have also performed tests with other, more complex models, like the mean-
field AdEx second order model [35, 34]. In this case, comparison can only be
done with TVB, since FAST_TVB only works with the reduced Wong Wang model and
does not allow new models to be easily incorporated. In both cases (i.e., TVB
and TVB C++), testing another model involves only the change of a single line
of code. In Figure 7(a) we show the excitatory firing rate of the 379-region
model over 10 seconds. Again, both results are numerically identical.

TVB | TVB C++
---|---
489.65 | 23.5
Table 2: Execution times for the mean-field AdEx first order model simulation
(in seconds).

Execution times are shown in Table 2. As can be seen, simulation times are
speeded up about 10x times with TVB C++. This takes into account that TVB
Python uses Numba to accelerate the core of the simulation and the computation
of the derivatives of the differential equations of the model. In TVB C++ we
have not applied any further specific optimization.

Thus, the last set of tests has been performed using the mean-field AdEx model
[38, 35, 34]. In this case, we have reproduced some of the results presented
in [34], specifically Figure 2 in the paper. The results are shown in Figure
7.

###  5.3 FRE Model

Montbrió and coworkers [36, 37], proposed a method to derive the Firing Rate
Equations (FREs) for networks of heterogeneous, all-to-all coupled quadratic
integrate-and-fire (QIF) neurons, which is exact in the thermodynamic limit,
i.e., for large numbers of neurons. Their formulation uses mean-field
equations to reveal a dual role of electrical synapses: First, they equalize
membrane potentials favoring the emergence of synchrony. Second, electrical
synapses act as “virtual chemical synapses”, which can be either excitatory or
inhibitory depending upon the spike shape. The original results offer a
precise mathematical explanation of the intricate effect of electrical
synapses in collective synchronization.

$\begin{gathered}\tau\frac{dr}{dt}=\frac{\Delta}{\tau\pi}+2ru-gr\\\
\tau\frac{du}{dt}=u^{2}+\bar{\eta}-(\pi\tau r)^{2}+(J+g\ln a)\tau
r\end{gathered}$ (2)

The parameters $\bar{\eta}$ and $\Delta$ are, respectively, the center and
width of the static distribution of inputs to the individual spiking neurons,
which is considered Lorentzian. Also, we have the recurrent inputs with their
synaptic weight $J$, and the electric synapses weighted by their respective
weight $g$.

This model can be extended to an Excitatory-Inhibitory (E-I) network of
quadratic integrate-and-fire neurons modeling a Pyramidal-Interneuronal
Network Gamma (PING) rhythm [54, 55], and then extend it again to a whole-
brain made of these E-I nodes, by formulating the system of equations:

$\begin{gathered}\tau^{e}\frac{dr_{j}^{e}}{dt}=\frac{\Delta^{e}}{\tau^{e}\pi}+2r_{j}^{e}u_{j}^{e}-g^{e}r_{j}^{e}\\\
\tau^{i}\frac{dr_{j}^{i}}{dt}=\frac{\Delta^{i}}{\tau^{i}\pi}+2r_{j}^{i}u_{j}^{i}-g^{i}r_{j}^{i}\\\
\tau^{e}\frac{du_{j}^{e}}{dt}={u_{j}^{e}}^{2}+\bar{\eta^{e}}-(\pi\tau^{e}r_{j}^{e})^{2}+I_{j}^{e}\\\
\tau^{i}\frac{du_{j}^{i}}{dt}={u_{j}^{i}}^{2}+\bar{\eta^{i}}-(\pi\tau^{i}r_{j}^{i})^{2}+I_{j}^{i}\\\
\tau^{N}\frac{dS_{j}^{ee}}{dt}=-S_{j}^{ee}+J^{N(ee)}r_{j}^{e}\\\
\tau^{N}\frac{dS_{j}^{ie}}{dt}=-S_{j}^{ie}+J^{N(ie)}r_{j}^{e}\\\
I_{j}^{e}=I_{ext}^{e}+\tau^{e}S_{j}^{ee}-J_{j}J^{G(ei)}\tau^{i}r_{j}^{i}+GJ^{A}\tau^{e}\sum_{k\neq
j}C_{kj}S_{k}^{ee}\\\
I_{j}^{i}=I_{ext}^{i}+\tau^{e}S_{j}^{ie}-J^{G(ii)}\tau^{i}r_{j}^{i}\end{gathered}$
(3)

Observe that, in this last expression, we have the long-range synaptic weight
$J^{A}$, identified with the AMPA current. Also, we introduced the slower NMDA
gating variable $S$, which modulates connections from the excitatory
population into both the excitatory and inhibitory populations, being
regulated by its respective time constant $\tau^{N}$. Finally, the inhibitory
connections have been associated with GABA, so we have labeled its constant as
$J^{G}(xi)$, with $x$ being the receiving population (excitatory $e$ or
inhibitory $i$). From this system, it is easy to recover the equations used by
Sanz Perl and co-workers [56] by removing the electrical coupling
($g^{e}=g^{i}=0$), removing the FIC mechanism, assuming no external
stimulation ($I_{ext}^{e}=I_{ext}^{i}$ = 0), considering that
$\tau^{e}=\tau^{i}=\tau^{m}$, equating $J^{long,e}=J^{ee}$, and considering
the steady state of all the gating variables $S^{e,i}$.

Once the full whole-brain model is implemented, to recover the original
formulation by Montbrió and coworkers [37], we need, first, to decouple all
populations by setting $G=0$, then decouple the inhibitory and excitatory
populations by choosing $J^{G}=0$, and finally approximate the behavior of the
gating variable $S$ by its steady state behavior. However, this last step
requires setting $\frac{dS}{dt}=0$, which is not possible without modifying
the code. However, depending on the integration algorithm used, we could set
$\tau^{N}=\Delta t$, which, for Euler-Maruyama will reduce to the correct
expression for $S$, and for other more advanced integrators (e.g., Runge-
Kutta) will lead to an average of the firing rate during the current time-
step. In particular, we have reproduced Figure 3 from the original paper,
using three values of the spike asymmetry parameter $a$: $0.25$, $1$, and $4$
(see Figure 8(a)).

Using the same settings, we can go beyond the results presented in the
original papers and extend the model to a whole-brain simulation. The results
for the behavior of one node, when coupled with all the others for values of
$G$ of $1/4$ and $1$, are shown in Figure 8. As expected, the behavior for a
given node changes significantly when coupled with the rest of the network,
showing significant oscillatory behavior.

###  5.4 BOLD simulations

Another basic feature of TVB C++ is the ability to use a BOLD monitor to
compute fMRI signals from the raw output of the models, such as the models
presented before. FAST_TVB also computes a BOLD signal from the reduced Wong
Wang model but is limited to the Balloon-Windkessel hemodynamic model. We have
compared BOLD signal generation in our model with both TVB 2.8.1 and FAST_TVB.
As already mentioned, since the details of the BOLD algorithm used are
different for both of these systems (see above), TVB C++ implements both of
them to perform fair comparisons.

####  5.4.1 BOLD monitor using HRF convolutions

TVB 2.8.1 implements, as its main option, a hemodynamic response function
(HRF) to compute the BOLD signal. In the current implementation, these HRF
kernels are implemented using a first-order Volterra kernel. To test BOLD
generation we have used a smaller parcellation to have shorter execution
times. In particular, we used the Desikan parcellation with 68 regions.
Execution times are shown in Table 3. In figure 5, middle row, we show the
plot of a BOLD signal over 2 minutes. As can be seen, both are numerically
identical.

TVB | TVB C++
---|---
250.895 | 4.41
Table 3: Execution times for the reduced Wang Wong BOLD generation, using the
Desikan parcellation (68 ROIs) and with the HRF convolution-based model.

From Table 3 we can observe a speedup of about 50x, which includes both the
Reduced Wong Wang model simulations, _and_ the BOLD simulation. This result is
mainly because the BOLD computations in TVB do not use any kind of
acceleration (e.g., using the _numba_ library) and are pure Python. In these
cases, the performance difference can be quite dramatic.

####  5.4.2 BOLD monitor using Balloon-Windkessel

We have also compared BOLD generation from electrical activity obtained from a
simulation with FAST_TVB. In this case, both FAST_TVB and TVB C++ use a
Balloon-Windkessel algorithm. We have used the same smaller parcellation to
have a meaningful comparison. Execution times are shown in Table 4. In Figure
5, first row, we show the plot of the BOLD signal over 2 minutes. Again, both
signals are numerically identical.

FAST_TVB | TVB C++
---|---
10.89 | 4.26
Table 4: Execution times for the reduced Wang Wong BOLD generation using the
Balloon-Windkesssel model, with the Desikan parcellation (68 ROI).

The 3x speedup found here could be explained because, with this parcellation
and this model, the BOLD computations take the bulk of the total time, and
these operations are vector and matrix intensive, and we strongly believe the
optimizations provided by the Eigen library are responsible.

Up to this point, all comparisons have been conducted within a single-threaded
process. However, several algorithms necessitate a fitting process of model
parameters, thereby requiring parallel execution of numerous simulations. In
this regard, both FAST_TVB and TVB C++ are equipped to leverage multi-core
CPUs for parallel execution. Nevertheless, for extensive parallel executions,
the utilization of GPUs proves to be more efficient.

The TVB project has introduced a GPU backend named RateML, which entails
rewriting the model derivative equations to enable their use within a CUDA
kernel. A comparison between RateML and TVB C++ is presented in Figure 9. The
findings demonstrate that for a small number of parallel executions (up to
180/190), TVB C++ exhibits superior performance, especially at lower numbers.
However, in scenarios requiring massive parameter sweeps, RateML surpasses TVB
C++ by an order of magnitude in terms of efficiency.

##  6 Summary and Discussion

As demonstrated, TVB C++ surpasses the original TVB Python implementation in
terms of speed and performance while preserving its flexibility. The primary
role of TVB C++ as an efficient Back-End for demanding computational tasks is
seamless and economical, thanks to its Python bindings tailored for ease of
use and adaptability. It is noteworthy that TVB C++ holds significant
potential for further optimization, leveraging available libraries such as
Intel MKL or CUDA. Additionally, unlike Python, which is hindered by the
Global Interpreter Lock, the C++ version of TVB C++ can effortlessly exploit
native threads. The current iteration of TVB C++ offers a straightforward
thread pool mechanism to capitalize on multicore architectures. TVB C++
constitutes a valuable addition to the TVB family owing to its high
compatibility, exceptional flexibility, and superior performance. It emerges
as an ideal backend for computationally intensive tasks, requiring no
additional expertise beyond TVB basics.

In conclusion, TVB C++ stands as a formidable asset in the realm of
computational neuroscience, offering a compelling blend of speed, flexibility,
and ease of use. Through meticulous optimization and seamless integration with
the renowned TVB platform, TVB C++ elevates the capacity for large-scale
simulations and heavy-duty computational tasks to unprecedented heights. Its
exceptional performance, outpacing the original TVB-Python implementation,
underscores its prowess as a backend solution for tackling complex
neuroscientific inquiries with efficiency and precision. Moreover, TVB C++’s
compatibility and intuitive Python bindings ensure a smooth transition for
users, requiring no additional expertise beyond fundamental TVB knowledge.

With the potential for further optimization and the ability to harness
advanced computing resources like supercomputers, TVB C++ emerges as an
indispensable tool for researchers seeking to push the boundaries of
computational neuroscience. By seamlessly marrying speed, flexibility, and
ease of use, TVB C++ paves the way for groundbreaking discoveries and
advancements in our understanding of the brain’s intricate dynamics.

## Code availability

All the code described in this paper can be found in a Github public
repository [52].

## Author contributions (CRediT)

PR, VJ and GD: conceptualization, methodology. IM, GZ, JF, MS and GP:
methodology. IM: development, software. IM, GP, GZ, GD, MS, PR: Visualization,
Investigation. PR, VJ, GP: Formal Analysis. IM, GZ, GP: Writing – original
draft. GP: Funding acquisition. All authors: Validation, Data curation,
Writing – review and editing

## Acknowledgements

We thank Michiel van der Vlag for assisting in the comparisons to Rate ML. We
thank Alain Destexhe, Damien Depannemaecker and Ernest Montbrió for their help
adapting the Mean-field AdEx and the ERP population models.

## Funding

IM and GP were partially funded by Grant PID2021-122136OB-C22 funded by
MCIN/AEI/10.13039/501100011033 and by ERDF A way of making Europe. This work
has been supported by the European Union’s Horizon 2020 research and
innovation program under Specific Grant Agreement No. 785907 (HBP SGA2) and
Specific Grant Agreement No. 945539 (Human Brain Project SGA3). This work was
supported by the Virtual Research Environment at the Charité Berlin – a node
of EBRAINS Health Data Cloud. Part of computation has been performed on the
HPC for Research cluster of the Berlin Institute of Health. PR acknowledges
support by EU Horizon Europe program Horizon EBRAINS2.0 (101147319), Virtual
Brain Twin (101137289), EBRAINS-PREP 101079717, AISN – 101057655, EBRAIN-
Health 101058516, Digital Europe TEF-Health 101100700; German Research
Foundation SFB 1436 (project ID 425899996); SFB 1315 (project ID 327654276);
SFB 936 (project ID 178316478; SFB-TRR 295 (project ID 424778381); SPP
Computational Connectomics RI 2073/6-1, RI 2073/10-2,RI 2073/9-1; DFG Clinical
Research Group BECAUSE-Y 504745852, PHRASE Horizon EIC grant; 101058240;
Berlin Institute of Health & Foundation Charité,

## Conflict of Interest Statement

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

## References

  * [1] Gewaltig, M.-O. & Diesmann, M.  Nest (neural simulation tool).  _Scholarpedia_ 2, 1430 (2007).
  * [2] Stimberg, M., Brette, R. & Goodman, D. F.  Brian 2, an intuitive and efficient neural simulator.  _eLife_ 8, e47314, DOI: 10.7554/eLife.47314 (2019).
  * [3] Dura-Bernal, S. _et al._ Netpyne, a tool for data-driven multiscale modeling of brain circuits.  _eLife_ 8, DOI: 10.7554/elife.44494 (2019).
  * [4] Hines, M. L. & Carnevale, N. T.  Neuron: A tool for neuroscientists.  _The Neuroscientist_ 7, 123–135, DOI: 10.1177/107385840100700207 (2001).
  * [5] Migliore, M., Cannia, C., Lytton, W. W., Markram, H. & Hines, M. L.  Parallel network simulations with neuron.  _Journal of Computational Neuroscience_ 21, 119–129, DOI: 10.1007/s10827-006-7949-5 (2006).
  * [6] Bekolay, T. _et al._ Nengo: a python tool for building large-scale functional brain models.  _Frontiers in Neuroinformatics_ 7, DOI: 10.3389/fninf.2013.00048 (2014).
  * [7] Dai, K. _et al._ Brain modeling toolkit: An open source software suite for multiscale modeling of brain circuits.  _PLOS Computational Biology_ 16, e1008386, DOI: 10.1371/journal.pcbi.1008386 (2020).
  * [8] Leon, P. S. _et al._ The virtual brain: a simulator of primate brain network dynamics.  _Frontiers in Neuroinformatics_ 7, DOI: 10.3389/fninf.2013.00010 (2013).
  * [9] Ritter, P., Schirner, M., McIntosh, A. R. & Jirsa, V. K.  The virtual brain integrates computational modeling and multimodal neuroimaging.  _Brain Connectivity_ 3, 121–145, DOI: 10.1089/brain.2012.0120 (2013).
  * [10] Sanz-Leon, P., Knock, S. A., Spiegler, A. & Jirsa, V. K.  Mathematical framework for large-scale brain network modeling in the virtual brain.  _NeuroImage_ 111, 385–430, DOI: 10.1016/j.neuroimage.2015.01.002 (2015).
  * [11] Schirner, M. _et al._ Brain simulation as a cloud service: The virtual brain on ebrains.  _NeuroImage_ 251, 118973, DOI: 10.1016/j.neuroimage.2022.118973 (2022).
  * [12] van der Vlag, M. _et al._ RateML: A code generation tool for brain network models.  _Frontiers in Network Physiology_ 2, DOI: 10.3389/fnetp.2022.826345 (2022).
  * [13] van der Vlag, M. _et al._ Vast tvb parameter space exploration: A modular framework for accelerating the multi-scale simulation of human brain dynamics, DOI: 10.48550/ARXIV.2311.13337 (2023).
  * [14] Guennebaud, G., Jacob, B. _et al._ Eigen v3.  http://eigen.tuxfamily.org (2010).
  * [15] Thomas Yeo, B. T. _et al._ The organization of the human cerebral cortex estimated by intrinsic functional connectivity.  _Journal of Neurophysiology_ 106, 1125–1165, DOI: 10.1152/jn.00338.2011 (2011).
  * [16] Allen, E. A. _et al._ Tracking Whole-Brain Connectivity Dynamics in the Resting State.  _Cerebral Cortex_ 24, 663–676, DOI: 10.1093/cercor/bhs352 (2012).  https://academic.oup.com/cercor/article-pdf/24/3/663/14099596/bhs352.pdf.
  * [17] Hansen, E. C., Battaglia, D., Spiegler, A., Deco, G. & Jirsa, V. K.  Functional connectivity dynamics: Modeling the switching behavior of the resting state.  _NeuroImage_ 105, 525–535, DOI: https://doi.org/10.1016/j.neuroimage.2014.11.001 (2015).
  * [18] Kringelbach, M. L., McIntosh, A. R., Ritter, P., Jirsa, V. K. & Deco, G.  The rediscovery of slowness: Exploring the timing of cognition.  _Trends in Cognitive Sciences_ 19, 616–628, DOI: 10.1016/j.tics.2015.07.011 (2015).
  * [19] Deco, G. _et al._ The arrow of time of brain signals in cognition: Potential intriguing role of parts of the default mode network.  _Network Neuroscience_ 7, 966–998, DOI: 10.1162/netn_a_00300 (2023).
  * [20] Jirsa, V. _et al._ The virtual epileptic patient: Individualized whole-brain models of epilepsy spread.  _NeuroImage_ 145, 377–388, DOI: 10.1016/j.neuroimage.2016.04.049 (2017).
  * [21] Makhalova, J. _et al._ Virtual epileptic patient brain modeling: Relationships with seizure onset and surgical outcome.  _Epilepsia_ 63, 1942–1955, DOI: 10.1111/epi.17310 (2022).
  * [22] Adhikari, M. H. _et al._ Decreased integration and information capacity in stroke measured by whole brain models of resting state activity.  _Brain_ 140, 1068–1085, DOI: 10.1093/brain/awx021 (2017).
  * [23] Rocha, R. P. _et al._ Recovery of neural dynamics criticality in personalized whole-brain models of stroke.  _Nature Communications_ 13, DOI: 10.1038/s41467-022-30892-6 (2022).
  * [24] Idesis, S. _et al._ A low dimensional embedding of brain dynamics enhances diagnostic accuracy and behavioral prediction in stroke.  _Scientific Reports_ 13, DOI: 10.1038/s41598-023-42533-z (2023).
  * [25] Stefanovski, L. _et al._ Linking molecular pathways and large-scale computational modeling to assess candidate disease mechanisms and pharmacodynamics in alzheimer’s disease.  _Frontiers in Computational Neuroscience_ 13, 54, DOI: 10.3389/fncom.2019.00054 (2019).
  * [26] Patow, G., Stefanovski, L., Ritter, P., Deco, G. & Kobeleva, X.  Whole-brain modeling of the differential influences of amyloid-beta and tau in alzheimer’s disease.  _Alzheimer’s Research & Therapy_ 15, DOI: 10.1186/s13195-023-01349-9 (2023).
  * [27] Gast, R. _et al._ Pyrates—a python framework for rate-based neural simulations.  _PLOS ONE_ 14, e0225900, DOI: 10.1371/journal.pone.0225900 (2019).
  * [28] Cakan, C., Jajcay, N. & Obermayer, K.  neurolib: A simulation framework for whole-brain neural mass modeling.  _Cognitive Computation_ 15, 1132–1152, DOI: 10.1007/s12559-021-09931-9 (2021).
  * [29] Schirner, M., Rothmeier, S., Jirsa, V. K., McIntosh, A. R. & Ritter, P.  An automated pipeline for constructing personalized virtual brains from multimodal neuroimaging data.  _NeuroImage_ 117, 343–357, DOI: 10.1016/j.neuroimage.2015.03.055 (2015).
  * [30] Schirner, M., McIntosh, A. R., Jirsa, V., Deco, G. & Ritter, P.  Inferring multi-scale neural mechanisms with brain network modelling.  _eLife_ 7, DOI: 10.7554/elife.28927 (2018).
  * [31] Vella, M. _et al._ libneuroml and pylems: using python to combine procedural and declarative modeling approaches in computational neuroscience.  _Frontiers in Neuroinformatics_ 8, DOI: 10.3389/fninf.2014.00038 (2014).
  * [32] Deco, G. _et al._ How local excitation–inhibition ratio impacts the whole brain dynamics.  _Journal of Neuroscience_ 34, 7886–7898, DOI: 10.1523/JNEUROSCI.5068-13.2014 (2014).  http://www.jneurosci.org/content/34/23/7886.full.pdf.
  * [33] Aquilué-Llorens, D., Goldman, J. S. & Destexhe, A.  High-density exploration of activity states in a multi-area brain model.  _Neuroinformatics_ DOI: 10.1007/s12021-023-09647-1 (2023).
  * [34] di Volo, M., Romagnoni, A., Capone, C. & Destexhe, A.  Biologically realistic mean-field models of conductance-based networks of spiking neurons with adaptation.  _Neural Computation_ 31, 653–680, DOI: 10.1162/neco_a_01173 (2019).
  * [35] Zerlaut, Y., Chemla, S., Chavane, F. & Destexhe, A.  Modeling mesoscopic cortical dynamics using a mean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons.  _Journal of computational neuroscience_ 44, 45—61, DOI: 10.1007/s10827-017-0668-2 (2018).
  * [36] Montbrió, E., Pazó, D. & Roxin, A.  Macroscopic description for networks of spiking neurons.  _Physical Review X_ 5, DOI: 10.1103/physrevx.5.021028 (2015).
  * [37] Montbrió, E. & Pazó, D.  Exact mean-field theory explains the dual role of electrical synapses in collective synchronization.  _Physical Review Letters_ 125, DOI: 10.1103/physrevlett.125.248101 (2020).
  * [38] El Boustani, S. & Destexhe, A.  A Master Equation Formalism for Macroscopic Modeling of Asynchronous Irregular Activity States.  _Neural Computation_ 21, 46–100, DOI: 10.1162/neco.2009.02-08-710 (2009).  https://direct.mit.edu/neco/article-pdf/21/1/46/818085/neco.2009.02-08-710.pdf.
  * [39] Carlu, M. _et al._ A mean-field approach to the dynamics of networks of complex neurons, from nonlinear integrate-and-fire to hodgkin–huxley models.  _Journal of Neurophysiology_ 123, 1042–1051, DOI: 10.1152/jn.00399.2019 (2020).
  * [40] Goldman, J. S. _et al._ A comprehensive neural simulation of slow-wave sleep and highly responsive wakefulness dynamics.  _Frontiers in Computational Neuroscience_ 16, DOI: 10.3389/fncom.2022.1058957 (2023).
  * [41] Buxton, R. B. & Frank, L. R.  A model for the coupling between cerebral blood flow and oxygen metabolism during neural stimulation.  _J. Cereb. Blood Flow Metab._ 17, 64–72 (1997).
  * [42] Friston, K. J., Mechelli, A., Turner, R. & Price, C. J.  Nonlinear responses in fMRI: the balloon model, volterra kernels, and other hemodynamics.  _Neuroimage_ 12, 466–477 (2000).
  * [43] Boynton, G. M., Engel, S. A., Glover, G. H. & Heeger, D. J.  Linear systems analysis of functional magnetic resonance imaging in human V1.  _J. Neurosci._ 16, 4207–4221 (1996).
  * [44] Polonsky, A., Blake, R., Braun, J. & Heeger, D. J.  Neuronal activity in human primary visual cortex correlates with perception during binocular rivalry.  _Nat. Neurosci._ 3, 1153–1159 (2000).
  * [45] Glover, G. H.  Deconvolution of impulse response in event-related BOLD fMRI1.  _Neuroimage_ 9, 416–429 (1999).
  * [46] Friston, K. J., Josephs, O., Rees, G. & Turner, R.  Nonlinear event-related responses in fmri.  _Magnetic Resonance in Medicine_ 39, 41–52, DOI: 10.1002/mrm.1910390109 (1998).
  * [47] Obata, T. _et al._ Discrepancies between BOLD and flow dynamics in primary and supplementary motor areas: application of the balloon model to the interpretation of BOLD transients.  _NeuroImage_ 21, 144–153, DOI: 10.1016/j.neuroimage.2003.08.040 (2004).
  * [48] Stephan, K. E., Weiskopf, N., Drysdale, P. M., Robinson, P. A. & Friston, K. J.  Comparing hemodynamic models with DCM.  _NeuroImage_ 38, 387–401, DOI: 10.1016/j.neuroimage.2007.07.040 (2007).
  * [49] Stephan, K. E. _et al._ Nonlinear dynamic causal models for fMRI.  _NeuroImage_ 42, 649–662, DOI: 10.1016/j.neuroimage.2008.04.262 (2008).
  * [50] Deco, G. _et al._ Whole-brain multimodal neuroimaging model using serotonin receptor maps explains non-linear functional effects of lsd.  _Current Biology_ 28, 3065–3074.e6, DOI: https://doi.org/10.1016/j.cub.2018.07.083 (2018).
  * [51] Deco, G. & Kringelbach, M. L.  Metastability and coherence: Extending the communication through coherence hypothesis using a whole-brain computational perspective.  _Trends in Neurosciences_ 39, 125–135, DOI: https://doi.org/10.1016/j.tins.2016.01.001 (2016).
  * [52] Martín, I.  TVB C++ Library.  https://github.com/neich/tvb-root-cpp (2022).
  * [53] Wong, K.-F. & Wang, X.-J.  A recurrent network mechanism of time integration in perceptual decisions.  _Journal of Neuroscience_ 26, 1314–1328, DOI: 10.1523/jneurosci.3733-05.2006 (2006).
  * [54] Dumont, G. & Gutkin, B.  Macroscopic phase resetting-curves determine oscillatory coherence and signal transfer in inter-coupled neural circuits.  _PLOS Computational Biology_ 15, e1007019, DOI: 10.1371/journal.pcbi.1007019 (2019).
  * [55] Reyner-Parra, D. & Huguet, G.  Phase-locking patterns underlying effective communication in exact firing rate models of neural networks.  _PLoS Comput. Biol._ 18, e1009342 (2022).
  * [56] Sanz Perl, Y. _et al._ The impact of regional heterogeneity in whole-brain dynamics in the presence of oscillations.  _Network Neuroscience_ 7, 632–660, DOI: 10.1162/netn_a_00299 (2023).

## Figure Labels

Figure 1: Overview of the original TVB and the TVB C++ system pipelines.
Figure 2: First level of detail of the architecture of the TVB system. Taken
from the TVB main documentation. Figure 3: Architecture of the TVB C++ system.
The architecture of TVB C++ revolves around two distinct interfaces tailored
to user interaction. The back-end consists of blocks that seamlessly interface
with various top application layers. The _data types_ serve as the lingua
franca between different components such as analyzers, visualizers,
simulators, and uploaders. These datatypes embody "active data," meaning that
when TVB C++ is configured with a database, data encapsulated data type
instances are automatically kept. Presently, the console interface operates
independently of the storage layer, storing results solely in memory.
Consequently, users employing the console interface must manually manage data
import and export operations. Figure 4: Main modules of TVB C++. On the left
branch are all the core simulation modules, i.e., model, coupling, integrator,
monitor, noise, BOLD, and the simulator itself. On the right, are the tools,
e.g., observables, and the numpy and scipy interfaces.

Figure 5: Different plots comparing different platforms with the Reduced Wong
Wang model. Top: Plot of a (direct) simulation with TVB, FAST_TVB, and TVB
C++. Middle: Plot of a BOLD signal with TVB C++ and TVB, using HRF
convolutions. Bottom, first two columns: Plot of a BOLD signal from a whole
brain simulation with TVB C++ and TVB, using the Balloon-Windkesssel model.
Observe that in all cases, all platforms compared in the same rows produce the
exact same result. Bottom right: Plot of the maximum firing rate of the neuron
population for the same model, with and without the FIC mechanism.

Figure 6: Different simulations comparing how adding noise changes the results
of the simulations with the Modified Balanced Excitation-Inhibition model.
Left column: 10s simulation comparing TVB and TVB C++. Middle column: 200s
simulation comparing TVB and TVB C++, both used to generate a BOLD signal.
Right column: two examples of a BOLD signal computed from the previous
simulation, using the same algorithm (standard TVB monitor convolution) with
both TVB and TVB C++. Observe that in all cases, despite the random number
generation is not exactly the same, all comparisons produce a very similar
global result. Noise used is $\sigma=1\mathrm{e}{-5}$ in both TVB and TVB C++.

(a)

(b)

(c)

Figure 7: Reproducing the results of mean-field AdEx module from di Volo et
al. [34]. (a) Plot of a whole brain simulation with TVB and TVB C++ for the
mean-field AdEx first-order model. (b) the stable firing rate of the model for
different values of the adaptation parameter $b_{e}$. (c) the firing rate
distribution is predicted by the model for $b_{e}=60$.

(a)

(b)

(c)

Figure 8: Using TVB C++to reproduce the result of a whole brain simulation
using the FRE model [37] (Figure 3 in the manuscript). (a) Three values of the
spike asymmetry parameter $a$ are used: $1/4$, $1$ and $4$. (b-c) Using a
Desikan parcellation with 70 ROI, setting the spike asymmetry parameter $a$ at
$1/4$ (b) and $1$ (c) and showing the variation in the firing rate when the
coupling parameter $G$ varies from $0$ to $5$.

(a)

(b)

Figure 9: On the left panel we see the execution times for RateML as the
number of parallel executions grows. From 1 to 2000 parallel runs, the time is
almost constant around 11.5 minutes. However, when the GPU limit is reached,
the time starts to behave linearly with the number of executions. On the right
panel, the execution times of RateML against TVB C++. We can observe a linear
behavior for TVB C++ almost since the beginning, being more efficient up to
180/190 parallel executions, when RateML outperforms it given the massively
parallel characteristic of modern GPUs.

⬇

1import tvbcpp

2import numpy as np

3import matplotlib.pyplot as plt

4import scipy.io as sio

5

6dt = 0.1 # Integration step in [ms]

7m = sio.loadmat("CNT_S01_structure.mat")["CNT_S01_structure"]

8m = m / np.max(m) * 0.2 # normallize matrix

9

10# Set the structural connectivity matrix

11tvbcpp.set_weights(m)

12n_rois = m.shape[0]

13

14print(f"Loaded connectivity matrix with {n_rois} number of ROIs")

15

16tvbcpp.set_model("Montbrio") # Choose the neuronal model

17

18# Configure the integrator.

19sigmas = np.array([0.0, 0, 0, 0, 1e-3, 1e-3])

20tvbcpp.set_integrator("EulerStochastic", dt, sigmas)

21

22# Configure the initial state of the model

23n_state_vars = 6

24initial_state = np.zeros((n_rois, n_state_vars))

25tvbcpp.set_initial_state(initial_state)

26

27tvbcpp.set_global_coupling(80.1) # Set the global coupling

28

29tvbcpp.set_model_parameter("delta_i", 1.18) # Homogeneous working point
[SanzPerl’23]

30

31T = 300000 # simulation time in [ms]

32tvbcpp.add_temporal_average_monitor(1.0, [0]) # temporal average monitor

33

34result = tvbcpp.run_sim(0, T) # Run the simulation

35

36data_raw = result[0] # get the result of the first monitor

37time_series = data_raw[1][:, :, 0]

38

39tr = 2000 # BOLD sampling rate in [ms]

40

41bold_2007b = tvbcpp.create_bold("Stephan2007b", tr) # Create a BOLD object.

42# Compute the BOLD signal from a time series of raw data

43ts_2007b, data_bold_2007b = bold_2007b.compute_bold(time_series, 1.0)

44

45# Subsample originial time series from Stephan2007b

46n_2007b = T // tr

47t_sub_2007b = np.zeros((n_2007b,))

48data_sub_2007b = np.zeros((n_2007b, data_bold_2007b.shape[1]))

49for n in range(n_2007b):

50 t_sub_2007b[n] = np.mean(ts_2007b[n*tr:(n+1)*tr])

51 data_sub_2007b[n, :] = np.mean(data_bold_2007b[n*tr:(n+1)*tr, :], axis=0)

Listing 1: Python bindings example. The output from this code gives the
results shown in Figure 8.

⬇

1# Plot for run_sim

2fig, axs = plt.subplots(nrows=3, ncols=1, figsize=(10,10), sharex=True)

3axs = axs.flatten()

4

5axs[0].plot(data_raw[0][5000:], data_raw[1][5000:,:,0])

6axs[1].plot(t_sub_2007b[2:], data_sub_2007b[2:, :])

7axs[2].plot(ts_tvb[2:], data_bold_tvb[2:, :])

8

9plt.show()

Listing 2: Code to plot the output from Listing 1


#  DESY-23-101 Repercussions of the Peccei-Quinn axion on QCD

G. Schierholz
Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany

###### Abstract

The axion, originally postulated by Peccei and Quinn to solve the strong CP
problem, has become of great interest in particle and astroparticle
phenomenology. Yet it has a problem. It is widely assumed that the axion
leaves the nonperturbative features of QCD, such as the axial anomaly and
chiral symmetry breaking, unscathed. This is, however, not the case. It turns
out that the anomalous coupling of the axion to the gauge bosons can be
integrated partially, leaving behind a path integral extending over
topologically trivial gauge potentials only. This has far reaching
consequences. We conclude that the Peccei-Quinn axion extension of the
Standard Model is not a viable theory.

##  1 Introduction

Quantum chromodynamics (QCD) decribes the strong interactions remarkably well,
from the shortest distances probed so far to macroscopic distances shaped by
the topological properties of the theory. Yet it faces a problem. The theory
allows for a CP-violating term $S_{\theta}$ in the action,

$S=S_{\rm QCD}+S_{\theta}\,.$ (1)

In Euclidean space-time it reads

$S_{\theta}=i\,\theta\,Q\,,\quad Q=\int d^{4}x\,P\,\in\,\mathbb{Z}\,,$ (2)

where $Q$ is the topological charge and

$P=-\frac{1}{32\pi^{2}}\,\epsilon_{\mu\nu\rho\sigma}\textrm{Tr}\left[F_{\mu\nu}F_{\rho\sigma}\right]\,,$
(3)

is the Pontryagin density. In this formulation $\theta$ enters as an arbitrary
phase, referred to as the vacuum angle, which can take values
$-\pi<\theta\leq+\pi$. Thus, there is the possibility of CP violation in the
strong interaction. A nonvanishing value of $\theta$ would result in an
electric dipole moment $d_{n}$ of the neutron. To date the most sensitive
measurements of $d_{n}$ are compatible with zero. The current upper bound is
$|d_{n}|<1.8\times 10^{-13}e\,\textrm{fm}$ [1], indicating that $\theta$ is
anomalously small at hadronic scales. This puzzle is referred to as the strong
CP problem, which is one of the most intriguing unsolved problems in the
elementary particles field.

A popular view is that the solution of the strong CP problem lies beyond QCD
and the Standard Model. The first instinct in such a situation is to propose a
new symmetry that suppresses, or rules out, CP-violating processes in the
strong interactions. Indeed, Peccei and Quinn [2] concocted a global chiral
U(1) symmetry, which effectively rotates the $\theta$ parameter into a
dynamical, CP-conserving field $\phi$, the axion. This is realized by
augmenting the $\theta$ term (2) with the axion action [3],

$S_{\theta}\rightarrow S_{\theta}+S_{\phi}=\int
d^{4}x\,\left[\frac{1}{2}\big{(}\partial_{\mu}\phi)\big{)}^{2}+i\left(\theta+\frac{\phi}{f}\right)\,P+L_{\rm
int}(\partial_{\mu}\phi\,J_{\mu}^{5})\right]\,,$ (4)

where $f$ is a scale factor, and the Lagrangian $L_{\rm int}$ contains the
derivative coupling to the light quark axial vector currents, collectively
denoted by $J_{\mu}^{5}$. We have put all CP violation in $S_{\theta}$ and
deal with real quark mass matrices only. To leading order in $1/f$ the
interaction with the $u$ and $d$ quarks is of the form

$L_{\rm
int}=\frac{f_{u}}{f}\,\partial_{\mu}\phi\,\bar{u}\gamma_{5}\gamma_{\mu}u+\frac{f_{d}}{f}\,\partial_{\mu}\phi\,\bar{d}\gamma_{5}\gamma_{\mu}d\,,$
(5)

where $f_{u}$ and $f_{d}$ are dimensionless coupling constants. For the most
general form of the interaction respecting the U(1) symmetry see, for example,
[4]. It is readily seen that the path integral is invariant under the global
change of integration variables $\phi(x)\rightarrow\phi(x)+\delta$, called
shift symmetry, which eliminates $\theta$ from the action. Beyond the original
idea, the axion has aroused great interest in particle and astroparticle
phenomenology. For example, the axion is considered a candidate for dark
matter. Axions also arise in fundamental theories such as string theory. This
has prompted laboratory, astrophysical and cosmological searches for axions
and axion-like particles [5].

In all this effort it has been assumed that the axion extension of the
Standard Model does not change the nonperturbative properties of QCD. That can
be doubted. There is strong evidence, both from purely theoretical
considerations and numerical investigations [6, 7, 8, 9], for a nontrivial
$\theta$ dependence of the theory, which is incompatible with the shift
symmetry. Already in 1979 ’t Hooft [7] predicted at least one phase transition
at $|\theta|<\pi$. This leads us to expect that the axion will seriously
affect the topological properties of the QCD vacuum, thus questioning the
consistency of the model. In the following we will treat the axion field as a
dynamical degree of freedom, whether it arises from the U(1) Peccei-Quinn
symmetry or from a more fundamental theory, and focus on QCD interactions.

##  2 Path integral

We may drop the $\theta$ term in (4) now. This leaves us with the action

$S=\int
d^{4}x\,\left[\frac{1}{2}\,\big{(}\partial_{\mu}\phi\big{)}^{2}+\frac{i}{f}\,\phi\,P+L(\partial_{\mu}\phi\,J_{\mu}^{5})\right]+S_{\rm
QCD}\,.$ (6)

The only known way to solve the theory directly is by numerically evaluating
the path integral on a finite volume $V$, for example on the lattice. The
generating functional for Green’s functions is

$\int\\!\mathcal{D}[A_{\mu}]\int\\!\mathcal{D}[\phi]\,\int\,[\cdots]\,\exp\Big{\\{}\\!-\\!\int
d^{4}x\,\Big{[}\frac{1}{2}\big{(}\partial_{\mu}\phi\big{)}^{2}+\frac{i}{f}\,\phi\,P+L(\partial_{\mu}\phi\,J_{\mu}^{5})\Big{]}-S_{\rm
QCD}\Big{\\}}\,,$ (7)

where $\,\int{\scriptstyle[\cdots]}$ stands for the integral over the fermion
fields, which we will not specify for the moment. Physical quantities, like
Green’s functions, will not depend on the mean value of $\phi$. It thus can be
integrated out. Writing

$\phi(x)=\hskip 2.84526pt{\bar{\hskip-2.84526pt\phi}}+\hskip
2.84526pt{\tilde{\hskip-2.84526pt\phi}}(x)\,,\quad\hskip
2.84526pt{\bar{\hskip-2.84526pt\phi}}\equiv\langle\phi\rangle=\frac{1}{V}\int
d^{4}x\,\phi(x)\,,$ (8)

we may divide the path integral over the axion field $\phi(x)$ into an
integral over $\hskip 2.84526pt{\bar{\hskip-2.84526pt\phi}}$ and a path
integral over the field $\hskip 2.84526pt{\tilde{\hskip-2.84526pt\phi}}(x)$
orthogonal to $\hskip 2.84526pt{\bar{\hskip-2.84526pt\phi}}$,

$\int\mathcal{D}[\phi]\,\cdots\,=\int_{\langle\hskip
2.84526pt{\tilde{\hskip-2.84526pt\phi}}\rangle=0}\\!\\!\\!\mathcal{D}[\hskip
2.84526pt{\tilde{\hskip-2.84526pt\phi}}]\int d\hskip
2.84526pt{\bar{\hskip-2.84526pt\phi}}\,\cdots\,.$ (9)

Under (8) the anomalous contribution to the action divides into two pieces,

$\int d^{4}x\,\phi(x)\,P(x)=\hskip
2.84526pt{\bar{\hskip-2.84526pt\phi}}\,Q+\int d^{4}x\,\hskip
2.84526pt{\tilde{\hskip-2.84526pt\phi}}(x)\,P(x)\,,$ (10)

while the kinetic term and the derivative couplings, which do not depend on
$\hskip 2.84526pt{\bar{\hskip-2.84526pt\phi}}$, are not affected. We thus can
carry out the integral over $\hskip 2.84526pt{\bar{\hskip-2.84526pt\phi}}$.
Due to the periodicity of $\exp\\{-i\hskip
2.84526pt{\bar{\hskip-2.84526pt\phi}}Q\\}$, the integration can be restricted
to the finite interval $[-\pi,\pi]$. It results in a Kronecker delta function,
$2\pi f\delta_{Q\,0}$, which restricts the gauge fields to the sector of zero
topological charge, $Q=0$. Up to a factor, we end up with the generating
functional

$\int_{Q=0}\\!\mathcal{D}[A_{\mu}]\int_{\langle\phi\rangle=0}\\!\mathcal{D}[\phi]\,\int\,[\cdots]\,\exp\Big{\\{}\\!-\\!\int
d^{4}x\,\Big{[}\frac{1}{2}\big{(}\partial_{\mu}\phi\big{)}^{2}+\frac{i}{f}\,\phi\,P+L(\partial_{\mu}\phi\,J_{\mu}^{5})\Big{]}-S_{\rm
QCD}\Big{\\}}\,,$ (11)

where we have changed $\hskip 2.84526pt{\tilde{\hskip-2.84526pt\phi}}$ to
$\phi$. Restricting the integral to $\langle\phi\rangle=0$ is common practice
in the calculation of the constraint effective potential [10]. It should also
be noted that in the continuum the gauge fields split naturally into quantum
mechanically disconnected sectors of fixed topological charge. On the lattice
that has been demonstrated in [9]. Hence, the limitation to $Q=0$ is just a
matter of selecting the appropriate sector. The path integral over $\phi$ can,
in principle, be sampled by a white noise in a finite element method.

We can go one step further now and cast the remainder of the anomalous
coupling $\phi P$ into the Lagrangian $L_{\rm int}$ by partial integration.
The Pontryagin density is known to be a total divergence,
$P=\partial_{\mu}\Omega_{\mu}^{\scriptscriptstyle(0)}$, where

$\Omega_{\mu}^{\scriptscriptstyle(0)}=-\frac{1}{8\pi^{2}}\,\epsilon_{\mu\nu\rho\sigma}\,\text{Tr}\left[A_{\nu}\partial_{\rho}A_{\sigma}+\frac{2}{3}A_{\nu}A_{\rho}A_{\sigma}\right]$
(12)

is the Chern-Simons density or $0$-cochain [11, 12]. For purposes of the
topological analysis we compactify space-time. We consider two common
compactifications $M$, the 4D ball $B_{4}$ [13] and the 4D torus $T_{4}$ [14,
15, 16]. The ball is defined by
$B_{4}=\\{x\,|\,\sum_{\mu=1}^{4}x_{\mu}^{2}\leq R^{2}\\}$, whereas the torus
is defined by $T_{4}=\\{x\,|\,0\leq x_{\mu}\leq L_{\mu}\,;\,\mu=1,\cdots,4\\}$
with periods $L_{\mu}$. The radius $R$ will be taken large, so as to
extrapolate the results to the infinite volume. There are no restrictions on
$L_{\mu}$. We impose the condition that the gauge fields be regular
everywhere, except at the boundaries. This can be achieved (subsequently) on
any gauge field configuration being generated by the functional (11). By
partial integration we thus get

$\int_{M}d^{4}x\,\phi\,P=\int_{\partial
M}d^{3}\sigma_{\mu}\,\phi\,\Omega_{\mu}^{\scriptscriptstyle(0)}-\int_{M}d^{4}x\,\partial_{\mu}\phi\,\Omega_{\mu}^{\scriptscriptstyle(0)}\,,$
(13)

where $\partial M$ denotes the boundary of $M$. The boundary of $B_{4}$ is the
3-sphere $S_{3}$. To keep the action finite, $(\partial_{\mu}\phi)^{2}$ must
vanish on $S_{3}$ as $R$ goes to infinity, which requires that
$\phi=\text{const.}+O(1/R^{2})$. From
$Q=\int_{S_{3}}\\!d^{3}\sigma_{\mu}\,\Omega_{\mu}^{\scriptscriptstyle(0)}=0$
then follows that
$\int_{S_{3}}\\!d^{3}\sigma_{\mu}\,\phi\,\Omega_{\mu}^{\scriptscriptstyle(0)}$
vanishes as well, up to corrections of $O(1/R^{2})$, which disappear in the
infinite volume limit. On the torus, periodicity of the action implies
periodic boundary conditions for the axion and periodic boundary conditions up
to a gauge transformation $\varv_{\mu}=g(x_{\mu}=L_{\mu})\,g^{-1}(x_{\mu}=0)$
for the gauge fields. The matrix functions $\varv_{\mu}$ can be categorized
into homotopy classes labelled by the topological charge
$Q\in\pi_{3}(\text{SU(3)})=\mathbb{Z}$. For $Q=0$ the transformation
$\varv_{\mu}$ is homotopically trivial, which means that it can be
continuously deformed to $\varv_{\mu}=1$ [17]. The result is $\int_{\partial
T_{4}}\\!d^{3}\sigma_{\mu}\,\phi\,\Omega_{\mu}^{\scriptscriptstyle(0)}=0$,
just like for $M=B_{4}$ in the infinite volume limit. The second term on the
right-hand side of (13) can be absorbed in the Lagrangian,

$L_{\rm
int}-\frac{i}{f}\partial_{\mu}\phi\,\Omega_{\mu}^{\scriptscriptstyle(0)}\rightarrow
L_{\rm int}\,,$ (14)

replacing $J_{\mu}^{5}\rightarrow
J_{\mu}^{5}-(i/f)\,\Omega_{\mu}^{\scriptscriptstyle(0)}$. This leads us to our
final result for the generating functional

$\int_{Q=0}\\!\mathcal{D}[A_{\mu}]\int\\!\mathcal{D}[\phi]\,\int\,[\cdots]\,\exp\Big{\\{}\\!-\\!\int
d^{4}x\,\Big{[}\frac{1}{2}\big{(}\partial_{\mu}\phi\big{)}^{2}\\!+L(\partial_{\mu}\phi\,J_{\mu}^{5})\Big{]}-S_{\rm
QCD}\Big{\\}}\,,$ (15)

where we can relax the constraint $\langle\phi\rangle=0$ now. It should be
noted though that in the numerical analysis special measures have to be taken
to keep $\langle\phi\rangle$ from running to infinity. A similar problem
arises in noncompact lattice QED [18, 19].

The Chern-Simons density $\Omega_{\mu}^{(0)}$ is a gauge variant quantity, and
one might wonder if the action is gauge invariant. Under a gauge
transformation $A_{\mu}\rightarrow g^{-1}(A_{\mu}+\partial_{\mu})g$ [12]

$\Omega_{\mu}^{\scriptscriptstyle(0)}\rightarrow\Omega_{\mu}^{\scriptscriptstyle(0)}-\frac{1}{8\pi^{2}}\,\varepsilon_{\mu\nu\rho\sigma}\,\partial_{\nu}\text{Tr}\left[\partial_{\rho}g\,g^{-1}A_{\sigma}\right]-\frac{1}{24\pi^{2}}\,\varepsilon_{\mu\nu\rho\sigma}\text{Tr}\left[\partial_{\nu}g\,g^{-1}\partial_{\rho}g\,g^{-1}\partial_{\sigma}g\,g^{-1}\right]\,.$
(16)

The last two terms can be written as a total divergence,
$\Omega_{\mu}^{\scriptscriptstyle(0)}\rightarrow\Omega_{\mu}^{\scriptscriptstyle(0)}+\partial_{\nu}\Omega_{\mu\nu}^{\scriptscriptstyle(1)}$,
where $\Omega_{\mu\nu}^{\scriptscriptstyle(1)}$ is the $1$-cochain [12]. This
changes the action by the amount
$-\int_{M}d^{4}x\,\partial_{\mu}\phi\,\partial_{\nu}\Omega_{\mu\nu}^{\scriptscriptstyle(1)}$.
Partial integration gives

$-\int_{\partial
M}d^{3}\sigma_{\nu}\,\partial_{\mu}\phi\,\Omega_{\mu\nu}^{\scriptscriptstyle(1)}+\int_{M}d^{4}x\,\partial_{\mu}\partial_{\nu}\phi\,\Omega_{\mu\nu}^{\scriptscriptstyle(1)}\,.$
(17)

The first term vanishes for the same reason $\int_{\partial
M}d^{3}\sigma_{\mu}\,\phi\,\Omega_{\mu}^{\scriptscriptstyle(0)}$ vanishes, the
second term vanishes by symmetry,
$\Omega_{\mu\nu}^{\scriptscriptstyle(1)}=-\Omega_{\nu\mu}^{\scriptscriptstyle(1)}$.

##  3 Implications for QCD

The Atiyah-Singer index theorem relates the topological charge $Q$ to the
number of right-handed ($n_{+}$) and left-handed ($n_{-}$) fermionic zero
modes, $Q=n_{+}-n_{-}$, where $n_{-}$ ($n_{+}$) can be interpreted as the
number of instantons (anti-instantons) in the vacuum. In two space-time
dimensions it has been shown that $Q=n_{+}$ for $Q\geq 0$ and $Q=-n_{-}$ for
$Q\leq 0$ [20, 21]. In four space-time dimensions the anomalous Ward identity
implies that on a background of topological charge $Q$ the fermion determinant
has exactly $|Q|$ zero modes [22]. Thus $|Q|=n_{+}+n_{-}$, in accord with the
index theorem and the two-dimensional result. This has been confirmed on the
lattice. For the first time Chiu [23] showed that either $Q=n_{+}$ or
$Q=-n_{-}$. In a recent dedicated investigation [24], following up on [9], we
have verified that $|Q|=n_{+}+n_{-}$ in the SU(3) Yang-Mills theory for
charges up to $|Q|=O(10)$ using the gradient flow. Thus, we can assume that
$n_{+}=n_{-}=0$ in the Peccei-Quinn model (and its generalizations [4]),
leaving behind a vacuum without any instantons and anti-instantons.
Furthermore, at infrared scales $\mu\lesssim\Lambda_{\rm QCD}$ the action of
the individual charge sectors reduces to $S_{\rm QCD}^{Q}\simeq 8\pi^{2}|Q|$
[9], which excludes any significant long-distance effects for charge $Q=0$. We
also could have divided the manifold $M$ into a set of open submanifolds
$M_{i}\subset M$ with $\hskip
2.84526pt{\bar{\hskip-2.84526pt\phi}}_{i}=\langle\phi\rangle_{M_{i}}$, which
covers $M$ up to a set of measure zero, and integrated $\int_{M_{i}}d^{4}x\,P$
to zero on any submanifold $M_{i}$. We will show now that the restriction to
gauge field configurations with trivial topology turns QCD into a different
theory [22]. Our main focus will be on the U(1) problem and chiral symmetry
breaking, which both are directly related to the topological features of the
QCD vacuum.

The QCD action has a chiral U(1) symmetry for small $u$ and $d$ quark masses
that is not realized in the real world. The violation of this symmetry can be
attributed to the anomaly in the regularized theory, which has played an
important historical role in establishing QCD as the theory of the strong
interactions. In particular, it ensures that the $\eta^{\prime}$ meson is
massive. Instantons can be seen to give rise to interactions that explicitly
break the symmetry [25]. In the case of trivial topology, the transformation

$q\rightarrow\exp\\{i\delta\gamma_{5}\\}\,q\,,\quad\bar{q}\rightarrow\bar{q}\,\exp\\{i\delta\gamma_{5}\\}$
(18)

with $q=u,d$ leaves the path integral unchanged,

$\int_{Q=0}\\!\\!\\!\mathcal{D}[A_{\mu}]\int\\!\mathcal{D}[q]\mathcal{D}[\bar{q}]\rightarrow\int_{Q=0}\\!\\!\\!\mathcal{D}[A_{\mu}]\,\exp\\{i2\delta
Q\\}\int\\!\mathcal{D}[q]\mathcal{D}[\bar{q}]=\int_{Q=0}\\!\\!\\!\mathcal{D}[A_{\mu}]\int\\!\mathcal{D}[q]\mathcal{D}[\bar{q}]\,,$
(19)

which indicates that the anomaly has no physical effect, and that the
regularized theory is invariant under the U(1) transformations (18). This
thwarts the solution of the U(1) problem [25]. Instead, we would find the
$\eta^{\prime}$ mass to be of the order of the $K$ and $\eta$ mass,

$m_{\eta^{\prime}}^{2}\approx 2m_{K}^{2}-m_{\eta}^{2}\,.$ (20)

In addition, the light mesons would appear in parity doublets of approximately
equal mass, but no such parity doubling is observed.

In QCD the topology of gauge fields is related to the chiral condensate by the
anomalous chiral Ward identity. For two mass-degenerate quark flavors of mass
$m>0$ the relation is

$-m\langle\bar{q}q\rangle=\left.\frac{1}{2V}\,\frac{\partial\ln Z}{\partial\ln
m}\right|_{\theta=0}=-\left.\frac{2}{V}\,\frac{\partial^{2}\ln
Z}{\partial\theta^{2}}\right|_{\theta=0}=2\frac{\langle Q^{2}\rangle}{V}\,,$
(21)

where $Z$ is the partition function, and $m=m_{u}=m_{d}$ and
$\langle\bar{q}q\rangle=\langle\bar{u}u\rangle=\langle\bar{d}d\rangle$. From
(21) we read off the familiar result [22]

$\chi_{t}=\frac{1}{2}\,\Sigma\,m\,,$ (22)

where $\chi_{t}=\langle Q^{2}\rangle/V$ is the topological susceptibility and
$\Sigma=-\langle\bar{q}q\rangle$, customarily taken in the chiral limit. As
the axion projects onto background gauge fields of zero topological charge, we
have

$\Sigma=0$ (23)

in the finite and infinite volume. QCD is thus deprived of one of its most
distinctive features, spontaneous chiral symmetry breaking. This result is in
broad agreement with instanton based models of the QCD vacuum, in which the
chiral condensate is proportional to a power of the density of instantons and
anti-instantons [26, 27]. Phenomenological estimates of the axion mass and
axion couplings assume $\Sigma\sim\Lambda_{\rm QCD}^{3}$. See, for example,
[4, 28]. Hence, most of the estimates need to be revised. Very likely, the
Peccei-Quinn model does not confine quarks and gluons either.

We do not expect any long-range correlations, as one might expect if the
integration of the gauge fields is restricted to fixed topological charge
[22]. Long-range correlations are characteristic of spontaneous symmetry
breaking, due to the presence of massless particles. The situation is special
here. As has been stressed above, integrating out $\hskip
2.84526pt{\bar{\hskip-2.84526pt\phi}}$ eliminates instantons and anti-
instantons altogether.

##  4 Conclusions

We conclude that the axion extension of the Standard Model is not a viable
theory. It does not only eliminate any dependence on the vacuum angle
$\theta$, but also removes instantons from the vacuum, thus litterally
‘throwing out the baby with the bath’. This has far reaching consequences. We
mentioned the ‘failure’ of the Adler-Bell-Jackiw anomaly and the restoration
of chiral symmetry. The importance of instantons for the long-distance
properties of the theory and for hadron phenomenology has been documented in
countless work.

A lattice simulation of the model would be desirable. So far this is not
possible though, due to the complex nature of the action. In [29] we have
performed a calculation for imaginary couplings $1/f$, which renders the
action real, thus making it suitable for Monte Carlo simulations. However, in
this case the charge is not integrated to $Q=0$, but to
$Q\rightarrow\pm\infty$, which we did not know at the time.

## Acknowledgment

I like to thank Meinulf Göckeler for carefully reading the manuscript and for
constructive criticism.

## References

  * [1] C. Abel, et al. Phys. Rev. Lett. 124 (2020) 081803 [arXiv:2001.11966 [hep-ex]].
  * [2] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38 (1977) 1440.
  * [3] R. D. Peccei, Lect. Notes Phys. 741 (2008) 3 [arXiv:hep-ph/0607268 [hep-ph]].
  * [4] J. E. Kim and G. Carosi, Rev. Mod. Phys. 82 (2010) 557 [erratum: Rev. Mod. Phys. 91 (2019) 049902] [arXiv:0807.3125 [hep-ph]].
  * [5] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner and K. A. van Bibber, Ann. Rev. Nucl. Part. Sci. 65 (2015) 485 [arXiv:1602.00039 [hep-ex]].
  * [6] C. G. Callan, R. F. Dashen and D. J. Gross, Phys. Rev. D 20 (1979) 3279.
  * [7] G. ’t Hooft, Nucl. Phys. B 190 (1981) 455.
  * [8] J. L. Cardy and E. Rabinovici, Nucl. Phys. B 205 (1982) 1.
  * [9] Y. Nakamura and G. Schierholz, Nucl. Phys. B 986 (2023) 116063 [arXiv:2106.11369 [hep-ph]].
  * [10] L. O’Raifeartaigh, A. Wipf and H. Yoneyama, Nucl. Phys. B 271 (1986) 653.
  * [11] W. A. Bardeen, Nucl. Phys. B 75 (1974) 246.
  * [12] M. L. Laursen, G. Schierholz and U. J. Wiese, Commun. Math. Phys. 103 (1986) 693.
  * [13] R. Jackiw, in Gauge theories of the eighties, R. Raitio and J. Lindfors (eds.), Lecture Notes in Physics 181, Springer (1983); MIT-CTP-1000.
  * [14] M. Lüscher, Commun. Math. Phys. 85 (1982) 39.
  * [15] M. Göckeler, M. L. Laursen, G. Schierholz and U. J. Wiese, Commun. Math. Phys. 107 (1986) 467.
  * [16] A. S. Kronfeld, M. L. Laursen, G. Schierholz and U. J. Wiese, Nucl. Phys. B 292 (1987) 330.
  * [17] R. Jackiw and C. Rebbi, Phys. Rev. Lett. 37 (1976) 172.
  * [18] M. Göckeler, R. Horsley, E. Laermann, P. E. L. Rakow, G. Schierholz, R. Sommer and U. J. Wiese, Nucl. Phys. B 334 (1990) 527.
  * [19] R. Horsley, Y. Nakamura, H. Perlt, D. Pleiter, P. E. L. Rakow, G. Schierholz, A. Schiller, R. Stokes, H. Stüben, R. D. Young and J. M. Zanotti, J. Phys. G 43 (2016) 10LT02 [arXiv:1508.06401 [hep-lat]].
  * [20] M. M. Ansourian, Phys. Lett. B 70 (1977) 301.
  * [21] N. K. Nielsen and B. Schroer, Nucl. Phys. B 127 (1977) 493.
  * [22] H. Leutwyler and A. V. Smilga, Phys. Rev. D 46 (1992) 5607.
  * [23] T. W. Chiu, Phys. Rev. D 58 (1998) 074511 [arXiv:hep-lat/9804016 [hep-lat]].
  * [24] Y. Nakamura, G. Schierholz and H. Stüben, in preparation.
  * [25] G. ’t Hooft, Phys. Rept. 142 (1986) 357.
  * [26] D. Diakonov, Proc. Int. Sch. Phys. Fermi 130 (1996) 397 [arXiv:hep-ph/9602375 [hep-ph]].
  * [27] T. Schäfer and E. V. Shuryak, Rev. Mod. Phys. 70 (1998) 323 [arXiv:hep-ph/9610451 [hep-ph]].
  * [28] G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro, JHEP 01 (2016) 034 [arXiv:1511.02867 [hep-ph]].
  * [29] Y. Nakamura and G. Schierholz, [arXiv:1802.09339 [hep-lat]].