# CommsVAE: Learning the brain’s macroscale communication dynamics using
coupled sequential VAEs

Eloy Geenjaar
TReNDS center
Georgia Institute of Technology
Atlanta, GA, USA
<EMAIL_ADDRESS>
&Noah Lewis
TReNDS Center
Georgia Institute of Technology
Atlanta, GA
&Amrit Kashyap
Department of Brain Simulation
Charite University Hospital
Berlin, Germany
&Robyn Miller
TReNDS Center
Georgia State University
Atlanta, GA
Vince Calhoun
TReNDS Center
Georgia Institute of Technology
Atlanta, GA

###### Abstract

Communication within or between complex systems is commonplace in the natural
sciences and fields such as graph neural networks. The brain is a perfect
example of such a complex system, where communication between brain regions is
constantly being orchestrated. To analyze communication, the brain is often
split up into anatomical regions that each perform certain computations. These
regions must interact and communicate with each other to perform tasks and
support higher-level cognition. On a macroscale, these regions communicate
through signal propagation along the cortex and along white matter tracts over
longer distances. When and what types of signals are communicated over time is
an unsolved problem and is often studied using either functional or structural
data. In this paper, we propose a non-linear generative approach to
communication from functional data. We address three issues with common
connectivity approaches by explicitly modeling the directionality of
communication, finding communication at each timestep, and encouraging
sparsity. To evaluate our model, we simulate temporal data that has sparse
communication between nodes embedded in it and show that our model can uncover
the expected communication dynamics. Subsequently, we apply our model to
temporal neural data from multiple tasks and show that our approach models
communication that is more specific to each task. The specificity of our
method means it can have an impact on the understanding of psychiatric
disorders, which are believed to be related to highly specific communication
between brain regions compared to controls. In sum, we propose a general model
for dynamic communication learning on graphs, and show its applicability to a
subfield of the natural sciences, with potential widespread scientific impact.

##  1 Introduction

Characterizing macroscale communication between brain regions is a complex and
difficult problem, but is necessary to understand the connection between brain
activity and behavior. The effect of neural systems on each other, or
connectivity, has also been linked to psychiatric disorders, such as
schizophrenia Friston (2002). Hence, gaining a deeper understanding of the
dynamics underlying the communication between brain regions is both from the
perspective of understanding how our brain facilitates higher-order cognition
and also to provide insight into and consequently how psychiatric disorders
arise. Static functional network connectivity (sFNC) and dynamic functional
network connectivity (dFNC), computed respectively as the Pearson correlation
between regional activation timeseries over the full scan duration (sFNC) or
on shorter sliding windows (dFNC), are among the most widely reported measures
of connectivity between brain regions Hutchison et al. (2013). These
approaches calculate the correlation between the timeseries of each brain
region to find their coherence, either across the full timeseries (sFNC) or by
windowing the timeseries and calculating the correlation within each window
(dFNC). Although extensions have been proposed Hutchison et al. (2013), along
with more complex connectivity measures, such as wavelet coherence Yaesoubi et
al. (2015) and Granger causality Roebroeck et al. (2005); Seth et al. (2015),
Pearson correlation remains the most prevalent measure of brain network
connectivity. And even the less commonly employed metrics have issues stemming
from some combination of sensitivity to noise, linearity, symmetry, or coarse
timescales. Most importantly, these approaches do not directly model the
communication, but rather analyze it post-hoc. The pursuit of instantaneous
communication between brain regions Sporns et al. (2021) and generative
approaches to model communication can potentially lead to models that closely
resemble effective connectivity Avena-Koenigsberger et al. (2018). An
important advantage of generative models is that they allow us to move away
from post-hoc inference from context-naïve metrics toward simulating
macroscale brain communication. We propose the use of recurrent neural
networks as generative models of communication on both simulated data and
functional magnetic resonance imaging (fMRI) data to validate and demonstrate
the specificity of our model. Our method complements connectivity metrics,
since it does more than quantify the aggregation of communication Avena-
Koenigsberger et al. (2018), but directly simulates it, and can thus be
analyzed using those same connectivity metrics.

Creating an accurate generative model of the macroscale communication dynamics
in the brain is hard, due to its complexity. However, there are some general
design principles the brain follows Sterling & Laughlin (2015). Generally, the
brain tries to minimize its energy use, which is likely also true for
communication in the brain. Macroscale communication is an energy-intensive
process and involves white matter tracts, which are essentially highways that
connect spatially separate parts of the brain. Hence, the amount of
information and the number of times information is sent should generally be
limited. The bits needed to convey are limited by the brain using sparse
coding, which at a lower scale is how neurons encode and communicate
information Olshausen & Field (2004). Although mechanistically different, due
to metabolic and volume constraints, macroscale communication presumably
exhibits strategies similar to sparse coding to efficiently transfer
information between neural populations Bullmore & Sporns (2012); Sterling &
Laughlin (2015). To incorporate this inductive bias into our model, we
regularize the communication from one region in our model to another using a
KL-divergence term to a Laplace distribution. The communication itself is
modeled as a normal distribution, but by minimizing its divergence with a
Laplace distribution we encourage sparser temporal communication. Encouraging
sparse temporal communication implies that information is only sent when
necessary, and lower rates of communication lead to reductions in the brain’s
energy requirements.

The goal of our model is to get a better idea of communication dynamics in the
brain. To evaluate whether our model is equipped to find the underlying
dynamics of a known generative model, we train it on simulated data. After we
show that our model finds the correct generative model from the simulated
data, we apply our model to neuroimaging data with fairly well-established
neural pathways. Although the pathways are well-established, we do not know
the exact ground truth communication underlying fMRI data. We do know that
communication dynamics depend on the task it is trying to perform, even
outside of explicit windows in which these tasks are performed. Hence, to
evaluate our model on fMRI data, we train it on three different tasks and use
logistic regression to predict what task a small window of timesteps is from
using the communication our model predicts. We compare our model to dFNC and
show that our method is more sensitive to the task, especially for small
windows where no explicit task is performed. This has implications for both
resting-state fMRI and task fMRI, since our model may be able to provide a
deeper look into how tasks shape macroscale communications in the brain.
Understanding task responses allows us to more accurately understand their
correlates to behavior and how communication in the brain shapes certain
behavior. Since our results indicate that our model is more sensitive to
communications even at resting-state that are still related to a task. Being
sensitive even outside of external driving of fMRI activity is impactful for
resting-state fMRI and psychiatric disorder research, where we expect
communication dynamics to be different between controls and patients.
Furthermore, since we establish our method on simulated data and then extend
it to real data, we provide an initial framework for scientists that may want
to extend our method into their area of expertise.

##  2 Method

#### Generative communication models

A generative communication model should meet the following requirements: it
should not be prone to noise, model the directionality of the communication,
and allow for instantaneous communication dynamics. Regarding noise,
macroscale communication in the brain is often assessed using fMRI, which is
prone to oscillatory sources of noise such as respiration Lund et al. (2006)
among others. Connectivity metrics based on correlations are more prone to
oscillatory noise since they are global and thus affect the synchrony between
fMRI signals in different regions at the same time Birn (2012). There are
generally, however, many reasons why fMRI signals can exhibit similar
activation patterns in multiple regions that may or may not always be
explained by communication Birn (2012). In our model, we model the more global
signal of each region separately from the communication dynamics. The general
trend of the signal is modeled by some initial inputs, and the sparsity
regularization on the communication helps push our model towards solutions
that model fewer spurious correlations. Secondly, the symmetry of many
connectivity metrics has been mentioned previously, and we solve this issue by
having each region explicitly send directional communication to each other
region. Thirdly, connectivity metrics often estimate the connectivity over a
temporal interval. For example, a correlation between two brain regions is
calculated over several timesteps $T$. This leads to a trade-off between
temporal resolution, and accuracy of the estimated correlation coefficient. In
fact, the maximum temporal resolution is two timesteps, since it is impossible
to calculate the variance over less than two samples. Our method bypasses this
issue by explicitly sending communication at each timestep. Thus, the temporal
resolution of the communication we estimate is the same as the temporal
resolution of the input signal.

#### Formal definition of our model

The data structure we assume in this work is a simple graph
$\mathcal{G}=(V,E)$, with $N$ vertices ($V$), and $M$ edges ($E$), as shown in
Figure 1A. Generally, we can assume a fully-connected graph, since our
proposed model can learn to not send any communication over an edge, meaning
it essentially ’trims’ that edge. Each vertex $v^{i}$ has a corresponding
timeseries $\textbf{x}^{i}_{1...T}$ with $T$ timesteps, this timeseries is
modeled as the following dynamical system.

$\dot{\textbf{x}}^{i}=\textbf{F}^{i}(\textbf{x}^{i}(t),\textbf{u}^{i}(t))$ (1)

Figure 1: This figure shows the composition of a single node from a graph on
the left. The communication and reconstruction RNNs are completely independent
of each other, the communication RNN sends signals to other nodes, whereas the
reconstruction RNN reconstructs the original timeseries from outside
communication and an initial condition $z_{0}$ (right).

One recent model that has been used to model this dynamical system is LFADs
Pandarinath et al. (2018), where they model the external inputs
$\textbf{u}(t)$ as samples from a normal distribution parameterized by a
controller that has access to the hidden state of the encoder and decoder. We
want to incorporate a more explicit inductive bias in our model, however. In
our case, the external inputs for vertex $v^{i}$ are communications $c^{\to
i}_{1...T}$ and external events $\textbf{e}(t)$. We use the notation $\to i$
to mean all incoming communication from neighbors of vertex $v^{i}$, $i\to$ to
mean all outgoing communication from vertex $v^{i}$ to its neighbors, and
$j\to i$ to mean communication from vertex $v^{j}$ specifically to vertex
$v^{i}$. The external events can be known external perturbations of the
dynamical system, but are not modeled using a controller like in LFADs
Pandarinath et al. (2018). We model the non-linear function $\textbf{F}^{i}$
for each vertex separately, using a recurrent neural network (RNN), with
initial inputs $\textbf{z}^{i}_{0}$. We find the initial inputs to the
dynamical systems using an encoder, with the following probability
distribution $\log
p(\textbf{z}^{i}|\textbf{x}^{i}_{1...T})=\log\mathcal{N}\left(\textbf{z}^{i};\bm{\mu}^{i}_{0},\bm{\sigma}^{i}_{0}\right)$.
This distribution is parameterized using a bi-directional RNN, and we sample
the initial distributions from it $\textbf{z}_{0}\sim
p(\textbf{z}^{i}|\textbf{x}^{i}_{1...T})$. Together with the decoder
$p(\hat{x}^{i}_{1...T}|\textbf{z}^{i}_{0},\textbf{c}^{\to
i}_{1...T},\textbf{e}_{1...T})$, which is parameterized as an RNN as well, the
encoder-decoder pair forms a sequential variational autoencoder (S-VAE) Zhao
et al. (2017); Kingma & Welling (2013), except we have added two extra inputs,
namely the communication from other vertices, and external events. The S-VAE
is visualized in Figure 1B. Lastly, the communication from one vertex to
another is also modeled as a recurrent neural network, with $\textbf{c}^{i}$
as a random variable. $\log
p(\textbf{c}^{i}|\textbf{x}^{i}_{1...t})=\log\mathcal{N}\left(\textbf{c}^{i};\bm{\mu}^{i\to
j}_{t},\bm{\sigma}^{i\to j}_{t}\right)$. Each communication from vertex $i$ to
vertex $j$ is thus sampled from this distribution $\textbf{c}^{i\to j}\sim
p(\textbf{c}^{i}|\textbf{x}^{i}_{1...t})$. This allows us to regularize the
communication distribution between each vertex to be a zero-mean unit-variance
multi-dimensional Gaussian, as with variational autoencoders Kingma & Welling
(2013). The communication network is visualized in Figure 1C. Since each
communication network predicts communication to another brain region
independently, we also add a router network that works similarly to the forget
gate in a long short-term memory (LSTM) network. The controller receives all
the incoming communication $\textbf{c}^{\to i}_{t}$ at timestep $t$, these are
fed through a multi-layer perceptron (MLP), which has the same number of
outputs as the size of the communication vector $\textbf{c}^{\to i}_{t}$. The
output is transformed through a sigmoid to be between [0, 1] and multiplied
with the incoming communication vector. This allows the router to block
information from certain a certain vertex that is already captured by
information sent from other vertices. The controller is visually represented
in Figure 1A.

With simplified notation, we thus have the following equations for the S-VAE.

$\displaystyle\textbf{h}^{i}_{T}$
$\displaystyle=\text{RNN}^{i}_{\text{Enc}}(\textbf{x}^{i}_{1...T})$ (2)
$\displaystyle\bm{\mu}^{i}_{0}$
$\displaystyle=\textbf{W}^{i}_{\mu}\textbf{h}^{i}_{T},\qquad\bm{\sigma}^{i}_{0}=\textbf{W}^{i}_{\sigma}\textbf{h}^{i}_{T},\qquad\textbf{z}^{i}_{0}\sim\mathcal{N}\left(\bm{\mu}^{i}_{0},\bm{\sigma}^{i}_{0}\right)$
(3) $\displaystyle\hat{\textbf{x}}^{i}_{1...T}$
$\displaystyle=\text{RNN}^{i}_{\text{Dec}}(\textbf{z}^{i}_{0},\;\textbf{c}^{\to
i}_{1...T},\;\textbf{e}_{1...T})$ (4) Then, for the communication network we
obtain the following equations. $\displaystyle\textbf{h}^{i}_{t}$
$\displaystyle=\text{RNN}^{i}_{\text{Comm}}(\textbf{x}^{i}_{\leq t})$ (5)
$\displaystyle\bm{\mu}^{i\to}_{t}$
$\displaystyle=\text{MLP}^{i}_{\text{Comm},\mu}(\textbf{h}^{i}_{t}),\qquad\bm{\sigma}^{i\to}_{t}=\text{MLP}^{i}_{\text{Comm},\sigma}(\textbf{h}^{i}_{t}),\qquad\textbf{c}^{i\to}_{t}\sim\mathcal{N}\left(\bm{\mu}^{i\to}_{t},\bm{\sigma}^{i\to}_{t}\right)$
(6) $\displaystyle\textbf{c}^{i\to}$ $\displaystyle=\textbf{c}^{\to
i}\odot\sigma\left(\text{MLP}^{i}_{\text{Router}}\left(\textbf{c}^{\to
i}\right)\right)`$ (7)

Where $\sigma(\cdot)$ is the sigmoid function, and
$\text{MLP}^{i}_{\text{Comm},\mu}$ and $\text{MLP}^{i}_{\text{Comm},\sigma}$
share every layer, except the last one.intense

We can then adapt the evidence lower-bound (ELBO) Kingma & Welling (2013) as
follows, where
$q_{\phi}^{i}\left(\textbf{z}^{i}_{0},\textbf{c}^{i\to}_{1...T}\;|\;\textbf{x}^{i}_{1...T}\right)$
refers to the approximate posterior, and
$p_{\theta}^{i}\left(\textbf{x}^{i}_{1...T}\;|\;\textbf{z}^{i}_{0},\textbf{c}^{\to
i}_{1...T},\textbf{e}_{1...T}\right)$ is the generative network, in our case
the recurrent decoder.

$\displaystyle\mathcal{L}(\theta,\phi,\textbf{x}_{1...T})$
$\displaystyle:=\sum_{i=1}^{N}\frac{1}{N}\bigg{[}-\text{D}_{\text{KL}}\left(q_{\phi}(\textbf{z}^{i}_{0}\;|\;\textbf{x}^{i}_{1...T})\;||\;p(\textbf{z})\right)$
(8)
$\displaystyle-\lambda\text{D}_{\text{KL}}\left(q_{\phi}(\textbf{c}^{i\to}_{1...T}\;|\;\textbf{x}^{i}_{1...T})\;||\;p(\textbf{c})\right)$
(9)
$\displaystyle+\mathbf{E}_{q_{\phi}^{i}\left(\textbf{z}^{i}_{0},\textbf{c}^{i\to}_{1...T}|\textbf{x}^{i}_{1...T}\right)}\left[\log
p_{\theta}^{i}\left(\textbf{x}^{i}_{1...T}\;|\;\textbf{z}^{i}_{0},\textbf{c}^{\to
i}_{1...T},\textbf{e}_{1...T}\right)\right]\bigg{]}$ (10)

Where $p(\textbf{z})$ and $p(\textbf{c})$ are the priors for the initial
inputs and the communications and $\lambda$ is the sparsity multiplier,
respectively. In our work, we take a zero-mean unit-variance multi-dimensional
Gaussian distribution for the initial inputs and a zero-mean $0.1$ standard
deviation Laplacian for the communications. This ensures that the model sends
communication sparsely to other regions.

#### Intuition and synthetic data

To gain some intuition about our model, we explain how our model works on the
simulated data we use in this work. Throughout its training process, each
region’s communication model will try to find temporal patterns in its region
that can help reconstruct another region’s timeseries. We assume that if this
is the case, namely that part of the signal is explained by a pattern in
another region at a previous timestep, then communication could have occurred.
A simple example is a communication between two regions, where the timeseries
for each region is uniform noise. Let us assume there is a linearly decaying
pulse in the first region at a time $T1$ $\text{Uniform}(0,70)$, that decays
until $T2~{}T1+\text{Uniform}(30,50)$. Once the pulse has decayed in the first
region at time $T2$, a pulse starts in the second region that linearly decays
at time $T3~{}T2+\text{Uniform}(30,50)$. These are exactly the parameters we
use to create our synthetic data. Now, if we only observe one of the two
timeseries, the pulses will seem sampled from a uniform distribution to the
decoder of that region. However, if we have access to both regions at the same
time and enough training data, our model should be able to figure out that the
pulse in the first region structurally precedes the pulse in the second
region. In our synthetic data, we add a third region that models a sinusoid
and has a pulse after the second region’s pulse that linearly decays at time
$T4~{}T3+\text{Uniform}(30,50)$. This region also has a different underlying
signal, namely a randomly shifted sinusoid (shifted between
$~{}\text{Uniform}(-20,20)$) to show that our model can separate global signal
from communication. A visual example of our synthetic data is shown in Figure
2A. We expect that the initial inputs $z_{0}$ models the global signal so that
the communication our model predicts is not affected by the sinusoidal signal.
As mentioned previously, fMRI data also exhibits global signal that may affect
connectivity metrics, such as physiological noise. Hence, we can infer whether
our model is affected by the global signal in our simulation data to make it
more realistic when transferred to biological data. The communication we thus
expect is shown in Figure 2B. Namely, the communication model in region 1
should learn to send a signal to region 2 when the pulse in its region has
almost decayed, and region 2’s communication model should learn to send a
pulse to region 3 once its pulse has almost decayed. The expected
communications in our synthetic data are shown in Figure 2B. To evaluate
whether our model finds the expected communication, we center a normal
distribution at the start of the pulse for regions $2$ and $3$. We then
correlate the communication from region $1$ to region $2$ and region $2$ to
region $3$ with that normal distribution for multiple shifts between $-20$ and
$2$. We take the maximum correlation and subtract two times the L1 sparsity of
the other communications from it to obtain a measure of similarity, where
higher is better.

Figure 2: This figure shows the simulated input signals on the left and the
communication we would expect our model to find. The communication would
happen from R1 to R2 (middle right) right before T2, and from R2 to R3 (bottom
right) right before T3.

#### Biological data

The fMRI data we use in this work is obtained from the open-access, under data
usage terms, HCP-1200 dataset Van Essen et al. (2013) for all subjects with
cortical surface timeseries data. We use three tasks in our work, the Motor,
Emotion Processing, and Language Processing tasks. The data is registered
using multimodal surface registration (MSM) Robinson et al. (2014; 2018), and
surfaces are constructed using Freesurfer Glasser et al. (2013); Fischl
(2012). We use the Yeo-17 atlas Yeo et al. (2011) to separate the cortical
surface data into $17$ separate regions. The spatial information within each
region is averaged, as is commonly done in communication analyses of the
brain. Each subject’s timeseries is then independently band-pass filtered
independently ($0.01-0.15$Hz) and linearly detrended using Nilearn package
Abraham et al. (2014). Each task has a different length, so we cut off tasks
at $176$ timesteps, which is the smallest common size of the timeseries. The
motor task consists of a visual cue, after which a subject has to start
performing a motor task. The motor tasks a subject performs are related to the
left hand, right hand, left foot, right foot, and tongue. For the emotion
task, the subject sees a visual cue and then has to match two shapes shown at
the bottom of the screen to a shape at the top. This sub-task is referred to
as the neutral task, and the other task has the subject match one of two faces
at the bottom of the screen with a face shown at the top. The faces either
look fearful or angry. The language task consists of a sub-task where the
participant has to perform arithmetic and push a button for the right answer.
The story sub-task presents the subjects with a short story of 5-9 sentences.
The timing of each task is shown, after convolving it with SPM’s hemodynamic
response function Penny et al. (2011), in Appendix B. We train one model on
all three tasks such that the prediction of the task is not influenced by the
fact that the model was specifically trained to model communication for that
task.

#### Implementation

Our model is implemented in Pytorch Paszke et al. (2017) and trained on an
internal cluster using single NVIDIA V100 and NVIDIA 2080 GPUs, with a batch
size of $8$, the Adam optimizerKingma & Ba (2014), a 1E-4 weight decay, a
learning rate of 0.01, $0.1$ epsilon, and $0.9$, $0.999$ as betas. We reduce
the learning rate when it plateaus using a scheduler. At each plateau, we
reduce the learning rate by multiplying it with $0.95$, with a patience of $6$
epochs, and a minimum learning rate of $1E-5$. We train each model for $200$
epochs, for the simulated data across $8$ seeds: ($42$, $1337$, $9999$, $50$,
$100$, $500$, $1000$, $5000$). The model is trained on the first three seeds
for the biological data, this is due to computational limitations. Since we
use multiple sparsities, the number of runs quickly explodes, and each model
on biological data takes about $5-6$ hours to train. The hidden size we use
for the encoder, decoder, and communication model are all $16$, and the latent
dimensions for $z_{0}$ is $32$. For the synthetic data we use sparsities:
($0.001$, $0.025$, $0.005$, $0.0075$, $0.01$, $0.025$, $0.05$, $0.075$, $0.1$,
$0.25$, $0.5$, $0.75$, $1.0$, $2.5$). For the real data we use a smaller
number of sparsities due to computation limitations, but also because we have
narrowed the sparsity down at this point to: ($0.001$, $0.005$, $0.01$,
$0.025$, $0.05$, $0.1$, $1.0$, $2.0$, $3.0$). In order to not mix subjects
across the training, validation, and test set, we split them based on subject.
Hence, a subject occurs with every task in one of the sets it was assigned to
and in no other sets.

#### Experiments

To test whether our model is interesting as a complement to connectivity
metrics, we evaluate our model in three different ways. To get a better idea
of how much the controller improves our model, we perform each experiment with
both models. We refer to the model without the controller as CommsVAE and the
model with the controller as Meta-CommsVAE in the rest of the text. Firstly,
we evaluate our model on the aforementioned synthetic data. We calculate the
correspondence between the expected communications and the communication that
our model predicts on a test set and look at how the sparsity parameter
affects the average correspondence to the test set. Note that when we use the
communication in these experiments, we take the mean of the Gaussian
distribution, since it is the most probably point. Then, we train both our
models on an fMRI dataset with three tasks, each model is trained on all three
tasks conjointly not to induce any task-specific model noise that can affect
the downstream prediction. Since we do not have a ground truth for the
communication of the fMRI signal, we determine the best model by looking at
the reconstruction error on the validation set with respect to the sparsity,
see Appendix A. We use the model that performs the best on average over all
the seeds for a certain sparsity parameter and uses it for task prediction. We
use $4$ different window sizes [$5$, $10$, $20$, and $50$] to compare our
model to dFNC. This allows us to determine whether our model has better
temporal resolution since the quality of the correlation will get lower with a
smaller dNFC window. For each window, we take the absolute value of our
model’s communication (since it can be both negative and positive) and sum it
across each window. Then, we train a logistic regression classifier with L1
sparsity on the training and validation set using the temporally summed
communication to predict what task each subject in the test set’s window is
from. This results in a classification accuracy across time for the four
windows, and we perform the same routine for the dFNC. Finally, we look at the
absolute sum of the communication for each task. We average the absolute
temporal sum over the subjects in the test set and visualize a graph of the
communication that exhibit more than $25\%$ of the maximum sum for each task.
The reason we need to take the absolute sum is that the model can send both
positive and negative communication.

##  3 Results

#### Synthetic data

The results for the synthetic data are shown in Figure 3. Generally, we see a
trend where higher sparsity parameters increase the correspondence to our
expected communication on the test set, up to a point where increased sparsity
suppresses any communication. There is a peak around $0.01$ in Figure 3A for
both models, after which the correspondence with the expected communication
decreases. This exactly supports our hypothesis that encouraging sparse
communication helps. To show what the communication looks like in the two
models we test, we take the test set synthetic data, predict the communication
for each, and use the synthetic generation parameters to align the
communication timeseries. After aligning the communication timeseries, we take
the average over the test set and visualize what the communication looks like
in Figure 3B. Right before the pulse ends in the region $1$, it communicates
with region $2$ before its pulse begins. The communication from region $2$ to
region $3$ starts exactly when the pulse is almost linearly decayed, right
before the pulse starts in region $3$. The main difference between our Meta-
CommsVAE and CommsVAE model is that in Figure 3B, the CommsVAE’s
communications between region $1$ and region $2$ are closer to the actual
pulse. Generally, we do not see a trend in terms of a difference between the
two CommsVAE models, at least for the synthetic data. This is probably because
the task itself is sparse and hence being able to increase sparsity using the
controller is unnecessary. An example of the reconstruction of validation data
by our model is shown in Appendix C.

Figure 3: This figure shows the performance of our models for different
sparsity parameters (x-axis) on the left and a comparison of our model’s
communications to the ground truth event times (black) on the right. The
values at which events happen are Uniformly sampled, but we shift each
timeseries in the test set to align the communication and events and then
average the communication over the test set.

#### fMRI data

The results to test both the specificity of our generative model and how
instantaneous the communication is are shown in Figure 4. We produce these
results for the best model we find based on the reconstruction error on the
validation set, see Appendix A. The best model is the Meta-CommsVAE model with
a sparsity of $0.025$. For the next two results, we will use this model. Our
results with this model clearly show that our method outperforms dynamic
functional connectivity (dFNC) in terms of predicting what task each window is
from. This is true for every window, and every window size. Furthermore, our
method only has a small decline in classification accuracy for smaller window
sizes, which supports our hypothesis that our model can find instantaneous
communication. Our method already starts being task-specific early on in time,
meaning that even before any of the tasks start, the communication our model
finds in the brain is task-specific, which is likely to be true if we knew the
ground truth of communication. One reason that the model increases in task
accuracy is that the communication model is an RNN and only has access to
points previously, so it does not have enough information to start sending
information earlier on in the timeseries. After a few windows, the accuracy of
the task increases and is then relatively stable. This is different for the
dFNC, which fluctuates, likely based on when the tasks start, but our model is
specific to the task essentially the same across all of the windows in the
timeseries. Thus, our model is even task-specific when there is relatively
little task signal from either of the tasks in a certain window. Another
result we can gather from Figure 4 is that the lower-triangular part of the
matrix is worse at predicting the task than the upper triangular part of the
communication matrix for each window. The reason we include these results is
that the dFNC is symmetric and the lower and upper triangular parts of the
matrix are thus the same. To make the comparison fair, we compare separately
compare the triangular parts of our communication matrix with dFNC. Since our
method is directed, we can also take the whole matrix, which outperforms both
the lower and upper triangular parts of that matrix, as expected.

Figure 4: This figure shows four panels, each with a different window size to
compute both our model’s communication density and dFNC. The y-axis shows the
classification accuracy of the task given only that window on the test set.

Our last experiment aims to evaluate whether our model finds explicitly
directed communication that is different across the tasks. The hemodynamic
response model over time for each sub-task is shown in Appendix B. The
absolute sum over the timeseries for each task is shown in Figure 5. There are
a few interesting results that align with known scientific processes. For
example, in all models, the communication from the earlier visual areas (V-A)
is stronger towards the later visual areas (V-B). Furthermore, the
communication from early visual areas (V-A) to the somatomotor regions (SM-B)
is the strongest in the motor task. Since the motor sub-tasks are preceded by
a visual cue, it is expected that the early visual areas send information to
the somatomotor areas, since those have to execute the task. Another
interesting result is that one of the control regions (C-A) in the motor task
is an important node in sending information to other brain regions. Generally,
most regions only send to other regions, without receiving any communications
from those regions, meaning that our model finds explicitly directed
communication. This aligns with our hypothesis that generative models should
be directed. This is especially obvious for the emotion task, where
communication from the early visual areas (V-A) to other regions is the most
dominant. The emotion task also contains a visual cue. The visual
communication in the language task is less strong compared to the other
regions, and there are many other communication dynamics, but since the math
sub-task is partly visual, we do expect some of the visual communication that
we observe in the figure. One of the large differences between the motor and
emotion, and language task is that the limbic system is much more involved in
the language task. We would potentially expect some limbic communication in
the emotion task as well since one of the sub-tasks pertains to seeing
angry/fearful faces. The limbic system composes the hippocampus and amygdala,
which are involved in the storage of memories and regulating emotional
memories, respectively. The reason that the language task could evoke
communication from the limbic system is thus that the stories evoke some
emotional memories in the subjects. The stories are from Aesop’s fables and
may have modulated some emotional memories in the subjects. However,
inferences like these need to be studied more robustly in future work.

Figure 5: This figure shows the absolute temporal sum over the predicted
communication dynamics from our model. The colors are scaled within each task,
and the colorbar goes from $0$ to the maximum of the absolute sum for each
task. The reason we do not show any numbers is because they are dimensionless.

##  4 Discussion

This work proposes a generative model for macroscale communication in the
brain that we show predicts directed, instantaneous, and task-specific
communication dynamics. We compare our model on a synthetic task that tests
the requirements we set out for our model and find that our model can
accurately uncover the communication dynamics we incorporate into the
synthetic data. Then, we train our two models on three fMRI tasks and find
that reconstruction error decreases with higher sparsity parameters,
indicating that the sparse prior on brain communication dynamics fits our
data. The sparsity of macroscale brain communication needs to be studied
further but aligns with findings in computational neuroscience Olshausen &
Field (2004) and the design principles of the brain Sterling & Laughlin
(2015). We use a model with the lowest reconstruction error to predict what
task small time windows are from in a test set and compare the classification
accuracy with a commonly used metric, dynamic functional connectivity (dFNC).
Our model performs much better with smaller window sizes, supporting our claim
that our model predicts instantaneous communication. Then, to understand the
added benefit of directed models, we look at the absolute sum of the
communication for each task and find that almost all, especially strong,
connections are directed. The directions and the strength of the communication
also align with neuroscientific priors we have about how regions affect each
other. Thus, our results confirm what we set out to incorporate into our model
from the beginning.

We want to expand on this model by applying it to resting-state fMRI data, and
explore whether our model can be used to more deeply understand psychiatric
disorders. Furthermore, since this model is generative, we can start to probe
the model and interpret the response functions when we synthetically increase
communication from one region to another at certain time points. We can design
communication patterns to test how communication in different regions affects
the signal in others.

### Acknowledgments

Data were provided by the Human Connectome Project, WU-Minn Consortium
(Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657)
funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for
Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at
Washington University. This material is supported by the National Science
Foundation under Grant No. 2112455 and the National Institutes of Health grant
#R01MH123610. Eloy Geenjaar was supported by the Georgia Tech/Emory NIH/NIBIB
Training Program in Computational Neural-engineering (T32EB025816).

## References

  * Abraham et al. (2014) Alexandre Abraham, Fabian Pedregosa, Michael Eickenberg, Philippe Gervais, Andreas Mueller, Jean Kossaifi, Alexandre Gramfort, Bertrand Thirion, and Gaël Varoquaux.  Machine learning for neuroimaging with scikit-learn.  _Frontiers in neuroinformatics_ , 8:14, 2014.
  * Avena-Koenigsberger et al. (2018) Andrea Avena-Koenigsberger, Bratislav Misic, and Olaf Sporns.  Communication dynamics in complex brain networks.  _Nature reviews neuroscience_ , 19(1):17–33, 2018\.
  * Birn (2012) Rasmus M Birn.  The role of physiological noise in resting-state functional connectivity.  _Neuroimage_ , 62(2):864–870, 2012.
  * Bullmore & Sporns (2012) Ed Bullmore and Olaf Sporns.  The economy of brain network organization.  _Nature reviews neuroscience_ , 13(5):336–349, 2012.
  * Fischl (2012) Bruce Fischl.  Freesurfer.  _Neuroimage_ , 62(2):774–781, 2012.
  * Friston (2002) Karl J Friston.  Dysfunctional connectivity in schizophrenia.  _World Psychiatry_ , 1(2):66, 2002.
  * Glasser et al. (2013) Matthew F Glasser, Stamatios N Sotiropoulos, J Anthony Wilson, Timothy S Coalson, Bruce Fischl, Jesper L Andersson, Junqian Xu, Saad Jbabdi, Matthew Webster, Jonathan R Polimeni, et al.  The minimal preprocessing pipelines for the human connectome project.  _Neuroimage_ , 80:105–124, 2013.
  * Hutchison et al. (2013) R Matthew Hutchison, Thilo Womelsdorf, Elena A Allen, Peter A Bandettini, Vince D Calhoun, Maurizio Corbetta, Stefania Della Penna, Jeff H Duyn, Gary H Glover, Javier Gonzalez-Castillo, et al.  Dynamic functional connectivity: promise, issues, and interpretations.  _Neuroimage_ , 80:360–378, 2013.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kingma & Welling (2013) Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Lund et al. (2006) Torben E Lund, Kristoffer H Madsen, Karam Sidaros, Wen-Lin Luo, and Thomas E Nichols.  Non-white noise in fmri: does modelling have an impact?  _Neuroimage_ , 29(1):54–66, 2006.
  * Olshausen & Field (2004) Bruno A Olshausen and David J Field.  Sparse coding of sensory inputs.  _Current opinion in neurobiology_ , 14(4):481–487, 2004.
  * Pandarinath et al. (2018) Chethan Pandarinath, Daniel J O’Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D Stavisky, Jonathan C Kao, Eric M Trautmann, Matthew T Kaufman, Stephen I Ryu, Leigh R Hochberg, et al.  Inferring single-trial neural population dynamics using sequential auto-encoders.  _Nature methods_ , 15(10):805–815, 2018.
  * Paszke et al. (2017) Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in pytorch.  2017\.
  * Penny et al. (2011) William D Penny, Karl J Friston, John T Ashburner, Stefan J Kiebel, and Thomas E Nichols.  _Statistical parametric mapping: the analysis of functional brain images_.  Elsevier, 2011.
  * Robinson et al. (2014) Emma C Robinson, Saad Jbabdi, Matthew F Glasser, Jesper Andersson, Gregory C Burgess, Michael P Harms, Stephen M Smith, David C Van Essen, and Mark Jenkinson.  Msm: a new flexible framework for multimodal surface matching.  _Neuroimage_ , 100:414–426, 2014.
  * Robinson et al. (2018) Emma C Robinson, Kara Garcia, Matthew F Glasser, Zhengdao Chen, Timothy S Coalson, Antonios Makropoulos, Jelena Bozek, Robert Wright, Andreas Schuh, Matthew Webster, et al.  Multimodal surface matching with higher-order smoothness constraints.  _Neuroimage_ , 167:453–465, 2018.
  * Roebroeck et al. (2005) Alard Roebroeck, Elia Formisano, and Rainer Goebel.  Mapping directed influence over the brain using granger causality and fmri.  _Neuroimage_ , 25(1):230–242, 2005.
  * Seth et al. (2015) Anil K Seth, Adam B Barrett, and Lionel Barnett.  Granger causality analysis in neuroscience and neuroimaging.  _Journal of Neuroscience_ , 35(8):3293–3297, 2015\.
  * Sporns et al. (2021) Olaf Sporns, Joshua Faskowitz, Andreia Sofia Teixeira, Sarah A Cutts, and Richard F Betzel.  Dynamic expression of brain functional systems disclosed by fine-scale analysis of edge time series.  _Network Neuroscience_ , 5(2):405–433, 2021.
  * Sterling & Laughlin (2015) Peter Sterling and Simon Laughlin.  _Principles of neural design_.  MIT press, 2015.
  * Van Essen et al. (2013) David C Van Essen, Stephen M Smith, Deanna M Barch, Timothy EJ Behrens, Essa Yacoub, Kamil Ugurbil, Wu-Minn HCP Consortium, et al.  The wu-minn human connectome project: an overview.  _Neuroimage_ , 80:62–79, 2013.
  * Yaesoubi et al. (2015) Maziar Yaesoubi, Elena A Allen, Robyn L Miller, and Vince D Calhoun.  Dynamic coherence analysis of resting fmri data to jointly capture state-based phase, frequency, and time-domain information.  _Neuroimage_ , 120:133–142, 2015.
  * Yeo et al. (2011) BT Thomas Yeo, Fenna M Krienen, Jorge Sepulcre, Mert R Sabuncu, Danial Lashkari, Marisa Hollinshead, Joshua L Roffman, Jordan W Smoller, Lilla Zöllei, Jonathan R Polimeni, et al.  The organization of the human cerebral cortex estimated by intrinsic functional connectivity.  _Journal of neurophysiology_ , 2011.
  * Zhao et al. (2017) Shengjia Zhao, Jiaming Song, and Stefano Ermon.  Towards deeper understanding of variational autoencoding models.  _arXiv preprint arXiv:1702.08658_ , 2017.

##  Appendix A fMRI data reconstruction error across sparsity

Figure 6 shows that the Meta-CommsVAE outperforms the CommsVAE on all
sparsities, but also has a clear dip around a certain sparsity. Since we do
not know the ground truth of the communication in the fMRI dataset like we did
with the synthetic data, the improvement in reconstruction error for higher
sparsities shows promise for hyperparameter searches on other datasets. Note
the reconstruction error is a squared error, summed over time and averaged
over the number of regions.

Figure 6: This figure shows the average reconstruction of both models on the
fMRI dataset, with standard deviation computed across the seeds. The Meta-
CommsVAE clearly outperforms the CommsVAE for all sparsities, and performs the
best at a sparsity of $0.025$.

##  Appendix B Tasks hemodynamic response

Figure 7 shows the hemodynamic response for each of the tasks for the $176$
timesteps we use in our model. The legend indicates what sub-task the
hemodynamic response corresponds to.

Figure 7: This figure shows the average hemodynamic response of the test set.

##  Appendix C Example reconstruction of the synthetic data

Figure 8 shows an example of the reconstruction of the synthetic data. The red
line is the input signal, and the blue line is the reconstruction. The first
row is the first region, the second row the second region, and the last row
the third region. The columns are four random timeseries from the validation
set. This shows that our model can still model the global signal, even for a
shifted sinusoid, with almost perfect communication dynamics, see Figure 3,

Figure 8: This figure shows an example of the reconstruction on the synthetic
data for our CommsVAE model with sparsity $0.01$. The red line is the input
signal, the blue line the reconstructed signal, each row corresponds to a
different region, and each column to a different timeseries from the
validation set.


conclusion applies as well to singlets.):

$M^{n}_{A}(Q^{2})=\left(\frac{\alpha_{c}(Q^{2})}{\alpha_{c}(Q^{2}_{o})}\right)^{d_{n}}M^{n}_{A}(Q^{2}_{0})\
,$ (6.2)

where $d_{n}\equiv-B_{n}/2\pi b>0$. If QCD is correct, and if $Q^{2}$ is large
enough so leading order in perturbation theory suffices and $O(1/Q^{2})$
corrections are negligible, then $\ln M^{n}_{A}(Q^{2})$ must lie on a straight
line when plotted versus $\ln[\alpha_{c}(Q^{2}_{0})/\alpha_{c}(Q^{2})]$ and he
slope must be $-d_{n}$. Such a plot is shown schematically in Fig. 25 for two
different targets with baryon numbers $A$ and $A^{\prime}$.

Figure 25: Rescaling a single moment.

At fixed $Q^{2}$, the EMC effect appears as the observation that
$M^{n}_{A^{\prime}}(Q^{2})<M^{n}_{A}(Q^{2})$ for $A^{\prime}>A$. On the other
hand, it is clear that there is a value of $Q^{2}$, call it $Q^{\prime\,2}$,
such that $M^{n}_{A^{\prime}}(Q^{\prime\,2})=M^{n}_{A}(Q^{2})$. As illustrated
in Fig. 25, $Q^{\prime\,2}<Q^{2}$. The value of $Q^{\prime\,2}$, in principle,
depends on $A,A^{\prime}$ and $Q^{2}$ and on $n$. Let us define
$\xi^{n}_{AA^{\prime}}(Q^{2})$ (not to be confused with the coordinate
$\xi_{\lambda}$ or the modified scaling variable $\xi$ of $\S 5$) so

$M^{n}_{A^{\prime}}(Q^{2})=M^{n}_{A}(\xi^{n}_{AA^{\prime}}(Q^{2})Q^{2})\ .$
(6.3)

From its definition

$\xi^{n}_{AA^{\prime}}(Q^{2})>1\ \ {\rm for\ }A^{\prime}>A$ (6.4)

and

$\displaystyle\xi^{n}_{AA^{\prime\prime}}(Q^{2})$
$\displaystyle=\xi^{n}_{AA^{\prime}}(Q^{2})\xi^{n}_{A^{\prime}A^{\prime\prime}}(Q^{2})$
$\displaystyle\xi^{n}_{AA^{\prime}}(Q^{2})$
$\displaystyle=1/\xi^{n}_{A^{\prime}A}(Q^{2})\ .$ (6.5)

The $Q^{2}-$dependence of $\xi^{n}_{AA^{\prime}}(Q^{2})$ is determined
entirely by QCD, and is independent of $A$ and $A^{\prime}$. Consider Eq.
(6.2) first for $A^{\prime}$ at $Q^{2}_{0}$ and $Q^{2}$

$M^{n}_{A^{\prime}}(Q^{2})=\left(\frac{\alpha_{c}(Q^{2})}{\alpha_{c}(Q^{2}_{o})}\right)^{d_{n}}M^{n}_{A^{\prime}}(Q^{2}_{0})$
(6.6)

and then for $A$ at $\xi^{n}_{AA^{\prime}}(Q^{2}_{0})Q^{2}_{0}$ and
$\xi^{n}_{AA^{\prime}}(Q^{2})Q^{2}$:

$M^{n}_{A}(\xi^{n}_{AA^{\prime}}(Q^{2})Q^{2})=\left(\frac{\alpha_{c}(\xi^{n}_{AA^{\prime}}(Q^{2})Q^{2})}{\alpha_{c}(\xi^{n}_{AA^{\prime}}(Q^{2}_{0})Q^{2}_{0})}\right)^{d_{n}}M^{n}_{A}(\xi^{n}_{AA^{\prime}}(Q^{2}_{0})Q^{2}_{0})\
.$ (6.7)

Now use Eq. (6.3) to eliminate all reference to moments

$\frac{\alpha_{c}(\xi^{n}_{AA^{\prime}}(Q^{2})Q^{2})}{\alpha_{c}(\xi^{n}_{AA^{\prime}}(Q^{2}_{0})Q^{2}_{0})}=\frac{\alpha_{c}(Q^{2})}{\alpha_{c}(Q^{2}_{0})}\
.$ (6.8)

To lowest order in $\alpha_{c}$,

$\alpha_{c}(\xi Q^{2})=\alpha_{c}(Q^{2})/(1+\alpha_{c}(Q^{2})b\ln\xi)\ ,$
(6.9)

which, together with Eq. (6.8), implies

$\xi^{n}_{AA^{\prime}}(Q^{2})=[\xi^{n}_{AA^{\prime}}(Q^{2}_{0})]^{\alpha_{c}(Q^{2}_{0})/\alpha_{c}(Q^{2})}\
.$ (6.10)

The extension to next order in $\alpha_{c}$ is quoted in [58]. Suppose
$Q^{2}_{0}<Q^{2}$, then $\alpha_{c}(Q^{2}_{0})/\alpha_{c}(Q^{2})>1$ and

$1<\xi^{n}_{AA^{\prime}}(Q^{2}_{0})<\xi^{n}_{AA^{\prime}}(Q^{2})\ \ {\rm for\
\ }Q^{2}_{0}<Q^{2}\ .$ (6.11)

The implication of this result is that a small change of $Q^{2}$-scale at low
$Q^{2}$ gets magnified into a large change when observed at large $Q^{2}$.

Eqs. (6.4), (6.5) and (6.10) summarize the properties of
$\xi^{n}_{AA^{\prime}}(Q^{2})$ which can be determined from general
considerations alone. The surprise came when Close, Roberts and Ross used this
method to analyze the EMC data and found that to a good approximation
$\xi^{n}_{AA^{\prime}}(Q^{2})$ appears to be independent of $n$:
$\xi^{n}_{AA^{\prime}}(Q^{2})\Rightarrow\xi_{AA^{\prime}}(Q^{2})$, at least
for the values of $n$ sensitive to the $x$-range of the EMC data. Actually CRR
did not construct moments but made an equivalent observation about the
structure function itself. Namely, if $\xi^{n}_{AA^{\prime}}(Q^{2})$ is
independent of $n$ then the structure functions themselves as functions of $x$
are related by a universal scale change in $Q^{2}$:

$\overline{F}^{A^{\prime}}_{2}(x,Q^{2})=\overline{F}^{A}_{2}(x,\xi_{AA^{\prime}}(Q^{2})Q^{2})\
.$ (6.12)

Eq. (6.12) has become known as “rescaling”. CRR were led to it by the
observation we referred to in $\S 5$, that the EMC effect resembles QCD
evolution. In fact, the EMC data are not in complete agreement with Eq.
(6.12). The excess at low $x$ is more than can be produced by the amount of
evolution that is required to fit the depletion at large $x$. Their analysis,
with $\xi_{DFe}\cong 2$ at $Q^{2}\approx$ 20 GeV2 is shown in Fig. 26.

Figure 26: Rescaling the EMC iron data [55].

The newer SLAC data on iron and deuterium data have a smaller enhancement at
low-$x$ and a lower cross over point (where $\overline{F}^{\rm
Fe}_{2}/\overline{F}^{D}_{2}=1$), both of which improve the agreement with the
rescaling analysis [57, 58]. The reader might wish to look back to Fig. 11 to
see the present state of the rescaling fits.

Several comments and caveats are in order:

  * 1.

QCD is subtle: Changing the scale creates quark-antiquark pairs. If, after all
the discussion of $\S 5$, this still seems unreasonable, perhaps it would help
to remember that a similar thing happens in a Bogoliubov transformation. By
redefining the vacuum, annihilation and creation operators get mixed up with
one another and a state which originally contained only particles, appears
after the transformation, to contain both particles and antiparticles.

  * 2.

Rescaling predicts that the EMC effect should vanish at $x\approx 0.2$ where
QCD evolution vanishes (see Fig. 23). The EMC data cross unity at $x\approx
0.35$ in disagreement with this prediction. Once again, however, the SLAC data
look better: $\overline{F}^{A}_{2}/\overline{F}^{D}_{2}$is definitely below
unity for $x>0.3$. A careful test must await better data at small-$x$.

  * 3.

Rescaling cannot work for $x\rightarrow 1$, or equivalently for
$n\rightarrow\infty$. At $x\rightarrow 1$, the structure function of the
nucleon vanishes but that of a nucleus does not. Eq. (6.12) fails to reproduce
this behavior. The reason for this failure will become clear soon.

  * 4.

Rescaling does not mean that an iron nucleus can be “mapped” to a nucleon by a
universal change of scale. The complete description of a nucleus requires
matrix elements of all twists. The gross structure of the EMC effect only
involves a few operators of twist-two.

### 6.2 Rescaling and Quark Models

In $\S 5.4$ we learned to associate a mass scale $\mu^{2}_{0}$ with the quark
model description of a hadron. It is the scale at which, approximately, the
hadron consists of valence quarks alone. It is an intrinsic characteristic of
each hadron. The measured structure function of the nucleon is consistant with
this notion with $\mu^{2}_{0}\cong 0.75$ GeV2. Clearly, any nucleus related to
the nucleon by rescaling, as in Eq. (6.12), also admits a valence quark
description, but at a shifted mass scale $\mu_{A}$:

$\mu^{2}_{A}=\frac{\mu^{2}_{0}}{\xi_{NA}(\mu^{2}_{0})}\ .$ (6.13)

Since $\xi_{NA}>1$, $\mu^{2}_{A}$ is less than $\mu^{2}_{0}$. This can be
checked by the following argument: At any fixed $Q^{2}$, the nucleus appears
more highly evolved than the nucleon (more pairs, softer valence quarks). So
to reabsorb all pairs it is necessary to devolve the nucleus further than the
nucleon.

It has become necessary to express Eq. (6.13) as a ratio of length scales. Let
us define $\lambda_{A}=c/\mu_{A}$, $\lambda_{N}=c/\mu_{0}$, then

$\lambda_{A}/\lambda_{N}=\mu_{0}/\mu_{A}=\sqrt{\xi_{NA}(\mu^{2}_{0})}\ $
(6.14)

$\xi(\mu^{2}_{0})$ is typically much smaller than $\xi(Q^{2})$, because of Eq.
(6.10). For example, $\xi_{\rm Fe\,N}(20\ {\rm GeV^{2}})\cong 2.0$, but
$\xi_{\rm Fe\,N}(\mu^{2}_{0})\cong 1.33$ making $\lambda_{\rm
Fe}/\lambda_{N}\cong 1.15$. Thus, an increase of only $\sim$ 15% in the
intrinsic length scale of the twist-two matrix elements for iron compared to
the nucleon can give rise to the EMC effect.

If there is anything fundamental about rescaling, rather than being merely an
accident, the basic relation must be defined at the intrinsic to the target,
$\mu^{2}_{0}$, not at some arbitrary $Q^{2}$ [55]. At $\mu^{2}_{0}$, however,
the structure function contains large contributions from higher twists, so an
equation like Eq. (6.12) with $Q^{2}$ replaced by $\mu^{2}_{0}$ cannot be
written down. Instead, we must write relations between matrix elements of
twist-two operators using the formalism of $\S 5.2$. We define

$\braket{P}{O_{n,a}^{(\mu^{2})\mu_{1}..\mu_{n}}}{P}=P^{\mu_{1}}...P^{\mu_{n}}\Theta^{A}_{n,a}(\mu^{2})$
(6.15)

in analogy to Eq. (5.31). Then the operator equivalent of Eq. (6.12) is

$A^{n-2}\Theta^{A}_{n,a}(\mu^{2}_{A})=A^{\prime\,{n-2}}\Theta^{A^{\prime}}_{n,a}(\mu^{2}_{A^{\prime}})\
.$ (6.16)

This is the basic rescaling between targets. The powers of $A$ and
$A^{\prime}$ are only kinematic. It leads to Eq. (6.12) at large $Q^{2}$
provided $\xi_{AA^{\prime}}(Q^{2})$ is obtained from
$\mu_{A}/\mu_{A}^{\prime}$ via QCD evolution (Eq. (6.10) plus higher order
improvements) and provided the moments can be reliably evolved from
$\mu_{A}^{2}$ to $Q^{2}$ using perturbative QCD. This turns out to be an
important proviso. In [52], we found that only the moments with $n<8$ could be
reliably devolved from large $Q^{2}$ to a $\mu^{2}_{0}$ as small as 1 GeV2.
The problem is that higher order (in $\alpha_{c}$) corrections typically
contain $\ln n$ factors making perturbation theory worse for large $n$.
Suppose, then, that Eq. (6.16) were exact for all $n$. Nevertheless, at large
$Q^{2}$ only the moments with $n<8$ would be likely to show uniform scaling.
This explains why rescaling fails near $x=1$ (as noted earlier) since large
$n$ moments probe exclusively large $x$. In fact, one can estimate the $x$
values for which rescaling should be reliable when observed at large $Q^{2}$
by using the Mellin transform relation [58].

$n\sim\ln\alpha_{s}/\ln x$ (6.17)

which for $n=8$ and $\alpha_{s}=0.2$ gives $x\sim 0.8$ as the upper limit of
reliability.

Turning this argument around one can see that exact uniform rescaling for all
$x$ at large $Q^{2}$ would be very hard to understand, since it would imply a
complicated and non-uniform relation among twist-two matrix elements with
$n>8$ at $\mu^{2}_{0}$.

The task for someone trying to understand the EMC effect from the point of
rescaling, then, is two-fold. First, one must explain why rescaling should be
uniform at the intrinsic scale, i.e. why should
$\lambda_{A}/\lambda_{A^{\prime}}$ be independent of $n$? Second, one must
predict the $A$ dependence of the rescaling parameter, i.e. of
$\lambda_{A}/\lambda_{N}$. Close, Roberts, Ross and I have argued that
rescaling is a rather natural prediction of quark models with only one length
scale. An example is the MIT bag model with only $u$ and $d$ quarks, which may
be taken to be massless, so the only dimensionful parameter is the bag
constant $B$. Such models are very close in spirit to QCD itself in which the
only dimensionful parameter is $\Lambda$. In a model like this, quarks
carrying momentum $p$ confined within a radius $\lambda$ transform into quarks
carrying momentum $p^{\prime}=(\lambda/\lambda^{\prime})p$ when the
confinement scale is changed to $\lambda^{\prime}$ \- the dimensionless
quantity $p\lambda$ is constant. The intrinsic scale $\mu^{2}$ is then
proportional to $p^{2}$, there being no other scale in the problem. It is hard
to turn this heuristic argument into a proof of uniform rescaling. That would
probably require a consistent formulation of perturbative QCD (including
renormalization) in a bag model, something which exists only in fragments
[62]. On the other hand, as Llewellyn Smith has noted, it seems clear that
other models such as the non-relativistic quark model have no hope of giving
rescaling unless the quark masses are assumed (rather unnaturally) to scale
with the inverse confinement radius [10]. A satisfactory theoretical
understanding of rescaling will have to await a more powerful QCD-based theory
of confinement. The second task $-$ determining the $A$ dependence of
$\lambda_{A}/\lambda_{N}$ $-$ is more straightforward. It is the subject of
the next Section.

### 6.3 A Dependence

It hardly needs saying that the EMC effect derives from the proximity of
nucleons within the nucleus, and that it would vanish if one could arrange
that that the nucleus was very dilute. Fortunately, nature has given us one
very dilute nucleus $-$ the deuteron $-$ and Bodek and Simon [15] have shown
that $\overline{F}^{D}_{2}(x,Q^{2})/F^{N}_{2}(x,Q^{2})$ is very close to
unity. It seems reasonable to assume, therefore, that the magnitude of the EMC
effect should be proportional to the probability that nucleons approach each
other or overlap within the nucleus.

Close, Roberts, Ross and I [57, 58] defined the simplest measure of this
effect we could imagine: we defined an “overlapping” volume for a nucleus
which is the integral over the nucleus of the two body density $\rho_{A}({\bf
r_{1}},{\bf r_{2}})$ multiplied by an overlap factor $V_{0}(|{\bf
r_{1}-r_{2}}|)$ which we took to be the overlapping volume of two spheres of
radius $a$:

$\displaystyle V_{0}(d)$
$\displaystyle=1-\frac{3}{4}\left(\frac{d}{a}\right)+\frac{1}{16}\left(\frac{d}{a}\right)^{3}{\
}$ $\displaystyle d\leq 2a$ $\displaystyle=0{\ \ \ \ \ \ }$ $\displaystyle
d>2a\ .$ (6.18)

Thus, the overlapping volume per nucleon is

$V_{A}=(A-1)\int d^{3}{\bf r_{1}}d^{3}{\bf r_{2}}\rho_{A}({\bf r_{1}},{\bf
r_{2}})V_{0}(|{\bf r_{1}-r_{2}}|)\ .$ (6.19)

$\rho_{A}({\bf r_{1}},{\bf r_{2}})$ is normalized to $\int d^{3}{\bf
r_{1}}d^{3}{\bf r_{2}}\rho_{A}({\bf r_{1}},{\bf r_{2}})=1$. Saturation of the
nuclear density at large $A$ implies $\rho_{A}\sim 1/A^{2}$. With this
behavior of $\rho_{A}$ and the finite integral of $V_{0}$ it is easy to see
that $V_{A}$ saturates at large $A$, i.e. $\lim_{A\rightarrow\infty}V_{A}=$
constant. The choice of a geometrical form for $V_{0}(|{\bf r_{1}-r_{2}}|)$
was in fact quite arbitrary. Any function which goes to unity as ${\bf
r_{1}-r_{2}}\rightarrow 0$ and to zero when $|{\bf r_{1}-r_{2}}|>2R_{\rm
nucleon}$ and which respects the three dimensional geometry of the problem
would do. We calculated the overlapping volume for nucleons with $a$ chosen so
$a_{\rm RMS}=0.9$ fm ($a_{\rm RMS}=\sqrt{\frac{3}{5}}a$, so $a=1.16$ fm).
$\rho_{A}({\bf r_{1}},{\bf r_{2}})$ was written in terms of the single
particle density $\rho_{A}({\bf r})$ and a correlation function $f({\bf
r_{1}-r_{2}})$:

$\rho_{A}({\bf r_{1}-r_{2}})=\rho_{A}({\bf r_{1}})\rho_{A}({\bf r_{2}})f({\bf
r_{1}-r_{2}})\ .$ (6.20)

$\rho_{A}({\bf r})$ was taken from experimental measurements of nuclear charge
densities. $f({\bf r_{1}-r_{2}})$ should, in principle, depend on $A$ but
there is little or no direct information on it from experiment. We took
$f({\bf r_{1}-r_{2}})$ from models of nuclear matter, the most realistic
probably being one based on a Reid soft core potential [60] shown in Fig. 27.

Figure 27: Choice of the correlation function $f(r)$: (a) no correlation; (b)
Fermi gas correlation; (c) Reid soft core correlation.

We then assumed that the effective confinement size in a nucleus $A$ will be
intermediate between that for an isolated nucleon and some limiting value
$\lambda_{\rm tot}$ associated with two totally overlapping nucleons. We
assumed a linear interpolation in $V_{A}$ leading to

$\frac{\lambda_{A}}{\lambda_{N}}=1+V_{A}\left(\frac{\lambda_{\rm
tot}}{\lambda_{N}}-1\right)\ .$ (6.21)

$\lambda_{\rm tot}$ might be viewed as a parameter. Instead, we estimated its
value from the bag model, where a spherical baryon number two system must have
volume at least twice as large as a single nucleon (otherwise it would be
stable against decay into two nucleons), so $R_{6}\geq 2^{1/3}R_{3}$. This led
us to take $\lambda_{\rm tot}/\lambda_{N}=2^{1/3}$, and to the values of
$\lambda_{A}/\lambda_{N}$ given in Table 2.

Table 2: Values of the confinement size relative to that for the free nucleon
for a range of nuclei. The three values (a), (b) and (c) correspond to three
choices of the correlation function $f(r)$ shown in Fig. 27. $\xi_{NA}$ (20
GeV2) is shown for the Reid correlation function. The others are similar.
Nucleus |  | $\lambda_{A}/\lambda_{N}$ |  | $\xi_{A}(Q^{2}=20)$
---|---|---|---|---
| (a) | (b) | (c) |
2D | 1.018 | 1.015 | 1.015 | 1.07
3He | 1.047 | 1.042 | 1.040 | 1.20
4He | 1.092 | 1.082 | 1.079 | 1.43
6Li | 1.054 | 1.045 | 1.045 | 1.23
7Li | 1.075 | 1.064 | 1.063 | 1.33
9Be | 1.088 | 1.074 | 1.074 | 1.40
12C | 1.124 | 1.105 | 1.104 | 1.60
16O | 1.128 | 1.109 | 1.108 | 1.63
20Ne | 1.122 | 1.104 | 1.104 | 1.60
27Al | 1.165 | 1.140 | 1.140 | 1.89
32S | 1.157 | 1.134 | 1.134 | 1.84
40Ca | 1.161 | 1.137 | 1.137 | 1.86
48Ca | 1.196 | 1.166 | 1.166 | 2.14
56Fe | 1.180 | 1.153 | 1.154 | 2.02
63Cu | 1.181 | 1.154 | 1.154 | 2.02
107Ag | 1.198 | 1.168 | 1.169 | 2.17
118Sn | 1.205 | 1.175 | 1.176 | 2.24
197Au | 1.229 | 1.196 | 1.195 | 2.46
208Pb | 1.220 | 1.188 | 1.188 | 2.37

These translate into values of $\xi_{NA}$(20 GeV2) which are much larger, and
these in turn give predictions for the EMC effect in a variety of nuclei. The
predictions of the rescaling model are compared with the SLAC data in Fig. 11.
The agreement is fine for $x<0.7$ above which Fermi motion becomes important.
One thing which is not obvious from Fig. 11 is that the data correlate well
with idiosyncracies in the periodic table. Fig. 28 shows predictions for
several $x$ values and $Q^{2}$ = 4.98 GeV2 compared with SLAC data.

Figure 28: $A$ dependence of the EMC effect at fixed $x$. The data are squared
with error bars. The predictions of rescaling are solid dots [58].

The fluctuations in the rescaling predictions reflect variations in nuclear
densities, e.g. 4He is more tightly bound than 3He, and generally follows the
data. Some predictions for the future are shown in Fig. 29.

Figure 29: $A$ dependence at $x=0.6$, $Q^{2}\approx 10$ GeV2 predicted by
rescaling.

Several comments are appropriate before leaving the discussion of $A$
dependence:

  * 1.

Llewellyn Smith has argued that the calculation of the $A$ dependence we have
given is much more general than the rescaling model. Any scheme which fits the
EMC effect in iron and in which the effect is linear in $V_{A}$ will agree as
well with the SLAC data. This is qualitatively true but the assumption of
linearity is non-trivial. We argued that $\lambda_{A}/\lambda_{N}$ should be
linear in the overlapping volume. Had we instead (erroneously) assumed
$\xi_{NA}(Q^{2}=20\ GeV^{2})$ to be linear in $V_{A}$ the $A$ dependence would
have come out wrong.

  * 2.

The rescaling model is not related to convolution models nor does it ascribe
the EMC effect to any particular exotic component in the nuclear wavefunction
(e.g., six quark bags). The scale change might originate from a change in the
size of individual nucleons, from quark percolation between nucleons, or from
multiquark or meson admixtures in the nuclear wavefunction. Personally, I
suspect that to the extent that these notions can be well-defined they will
turn out to correspond to the same underlying physics — partial deconfinement.

  * 3.

The rescaling fit to the $A$-dependence depends on several parameters: the
nucleon radius $a$, the formula chosen to relate $V_{A}$ to
$\lambda_{A}/\lambda_{N}$, the scale $\mu_{0}$ and the choice of two nucleon
correlation function. In fact, none of these were treated as free parameters
in [58]. Instead, they were fixed at “reasonable” values from other aspects of
hadron dynamics. Of course, the precise choice of parameters is not
particularly central to the explanation of the EMC effect.

  * 4.

The value of $V_{A}$ is very large for large $A$. It is $\approx$ 0.72 for
208Pb with the Reid correlation function. It has been remarked that such a
large value of $V_{A}$ is “unreasonable” and in some sense contradicted by the
many successes of conventional nuclear physics. The overlapping volume of
nucleons with $a_{\rm rms}$ = 0.9 fm is very large. That is a fact, not a
shortcoming of the rescaling model. As for the idea that large $V_{A}$ is
inconsistent with nuclear physics in general - I believe it is an unwarranted
concern: It costs $\approx$ 1 GeV (the string tension) to separate colored
sources by l fm. It even costs 300 MeV to flip a quark spin. These energetic
considerations, not naive classical, geometrical considerations determine the
quantum mechanics of nuclei.

### 6.4 Shadowing

It has long been expected that at very low values of $x$ deep inelastic
electron scattering from nuclei would behave like hadron scattering from
nuclei and exhibit “shadowing” [61]. At high energies hadron nucleus total
cross sections grow like $A^{2/3}$. The simple explanation of this effect is
that the incoming hadron doesn’t “see” the nucleons at the back of the
nucleus. They are in the shadow of those in front. The cross section then
grows like $\pi R^{2}\sim A^{2/3}$. Real photon-nucleus interactions are
shadowed [62]. This finds a simple explanation in the framework of vector
meson dominance: the hadronic interactions of a real photon are well-
approximated by supposing the photon converts with probability
$\sim\\!\\!\alpha$ into a vector meson ($\rho$, $\omega$ or $\phi$) which then
interacts hadronically and experiences shadowing. Vector meson dominance fails
to explain inelastic lepton scattering at large $Q^{2}$, or at least an
infinite tower of vector states with precisely tuned couplings is required,
but the prejudice that shadowing occurs there too is strong.

Figure 30: Cartoons representing the $\xi^{3}$ values probed by deep inelastic
lepton scattering at low-$x$ from an iron target.

My own thoughts on shadowing are still in flux. Since much of the future
experimental work in this field will be carried out in the kinematic regime
where shadowing might be important, I would like to outline the problem here.
Perhaps some reader will solve it! The kinematic relation

$\xi^{3}<1/Mx$ (6.22)

tells us that for very small values of $x$ the struck quark may propagate over
very large distances in the target. For $x\approx 0.5$, $\xi^{3}\leq 0.4$ fm,
but for $x\approx 0.05$, $\xi^{3}\leq 4$ fm and for $x\approx 0.005$,
$\xi^{3}\leq 40$ fm! The parton model diagram dominant at large $Q^{2}$ is
shown relative to an iron nucleus in Fig. 30 where the struck quark is
associated with a single nucleon.

Other possibilities represent EMC-like corrections. The question is whether a
quark propagating over such distances is or is not strongly absorbed in the
nuclear medium. Hadron interactions at very high energies are predominantly
absorptive, but the state which is propagating in Fig. 30 is not an ordinary
hadron. In particular, the propagating quark, which is a color triplet is
always close to an antiquark [63] which is a color antitriplet that may
neutralize some or most of its strong interactions. Furthermore, the invariant
mass of the $q\bar{q}$ pair is $-Q^{2}$.

Suppose for the sake of discussion, we characterize quark propagation in the
nuclear medium (accompanied by an antiquark as in Fig. 30) by an absorption
length $\lambda$ which may depend on $Q^{2}$. For a given nucleus when
$1/Mx<\lambda(Q^{2})$ there is little shadowing. Shadowing sets in when
$1/Mx\sim\lambda(Q^{2})$, but when $1/Mx$ exceeds twice the radius of the
nucleus ($R_{A}$) shadowing saturates: decreasing $x$ further does not put
more matter along the propagating quark’s path. So there are two $x$-values
characterizing shadowing, $x_{0}(Q^{2})=1/M\lambda(Q^{2})$ marking its onset
and $x(A)=1/2MR_{A}$ marking its saturation. For $x<x(A)$, only the front skin
of the nucleus down to a depth $\lambda(Q^{2})$ participates in the
scattering, so shadowing increases with $A$ but goes away as
$\lambda(Q^{2})\rightarrow\infty$. All of this can be summarized by the purely
phenomenological formula

$\overline{F}^{A}_{2}(x,Q^{2})=\frac{F^{N}_{2}(x,Q^{2})}{\left[1+\frac{x_{0}(Q^{2})}{x+x(A)}\right]}$
(6.23)

in which for simplicity I have ignored all other nuclear effects such as
rescaling and Fermi motion. Eq. (6.23) does not satisfy the quark number sum
rule. It is not intended to be valid at all $x$, only for $x\sim 0$ where
shadowing may be important. According to Eq. (6.23) there is an asymptotic
shadowing curve for nuclear matter ($A\rightarrow\infty$, $x(A)\rightarrow
0$), shown for example in Fig. 31. Finite nuclei track along
$\overline{F}^{\infty}_{2}(x,Q^{2})$ until $x\sim z(A)$ below which they
depart above $\overline{F}^{\infty}_{2}$.

A determination of $\lambda(Q^{2})$ is crucial [64]. Unfortunately, it is
largely unknown. Naively, one might expect $\lambda(Q^{2})$ to be a typical
meson mean free path, corresponding to a cross section of order 30 mb. This is
too naive. The propagating $q\bar{q}$ pair in Fig. 30 must have invariant mass
$-Q^{2}$. This can be generated if both quark and antiquark have transverse
momentum of order $\frac{1}{2}\sqrt{Q^{2}}$ or by a large mismatch in their
longitudinal momentum. In the former case, the transverse dimension of the
$q\bar{q}$ system is very small, $\sim 1/\sqrt{Q^{2}}$. The system looks like
a small color dipole and has a small absorption cross section. In the latter
case, the transverse dimensions are not small, color screening is not so
effective and the absorption cross section may be large. The importance of
shadowing depends on which region of phase space dominates. This is not yet
known. But at least some of the shadowing seen at $Q^{2}=0$ should disappear
at large $Q^{2}$. The behavior of $\overline{F}^{A}_{2}(x,Q^{2})$ at small $x$
expected if shadowing indeed disappears at large $Q^{2}(x(Q^{2})\rightarrow
0$) as shown in Fig. 31.

Figure 31: Scenarios for shadowing. (a) and (b) compared to two different
values of $\lambda(Q^{2})$ with $\lambda_{a}>\lambda_{b}$.

The discussion in this section has been very qualitative and phenomenological.
Eq. (6.23) should be taken with a grain of salt. There is much to do on the
subject of shadowing and little of substance to report here. Nevertheless, I
thought it might be appropriate to end by whetting the reader’s appetite for
the next round of experiments.

## References

  * [1] M. B. Johnson and A. Picklesimer, eds., Relativistic Dynamics and Quark Nuclear Physics, (Wiley-Interscience, New York, 1986).
  * [2] J.J. Aubert, et al., Phys. Lett. 123B, 275 (1983).
  * [3] Throughout these lectures we use the metric $g_{00}=-g_{ii}=1$, covariantly normalized states, $\braket{p}{p^{\prime}}=(2\pi)^{3}2E\delta^{3}({\bf p}-{\bf p^{\prime}})$, and the $\gamma$-matrices and other conventions of J.D. Bjorken and S.D. Drell, Relativistic Quantum Mechanics (McGraw-Hill, New York, 1964).
  * [4] In fact, the operator product in Eq. (1.8) may be so singular as $\xi_{\lambda}\rightarrow 0$ that the time ordered product is not Lorentz-invariant. In this event the proper definition of $T_{\mu\nu}$ makes use of the $T^{\star}$-product which differs from the $T$-product by local terms (distributions at $\xi_{\lambda}$ = 0) which restore Lorentz invariance. It is possible to derive paradoxical (and incorrect) results by ignoring this sublety. For a discussion, see the book by Itzykson and Zuber referenced in the preface.
  * [5] J.D. Bjorken, Phys. Rev. 148, 1467 (1966); 179, 1547 (1969).
  * [6] For an introduction to the methods of dispersion theory, see G. Barton, Dispersion Techniques in Field Theory (W.A. Benjamin, New York, 1965).
  * [7] L.D. Landau, Nucl. Phys. 13, 181 (1959); J.D. Bjorken, Ph.D. dissertation, Stanford University, 1959.
  * [8] For introductory treatments of the $P_{\infty}$ approach to deep inelastic scattering see the books by Close and Feynman referenced in the preface.
  * [9] M.J. Lighthill, Introduction to Fourier Analysis and Generalized Functions (Cambridge University Press, Cambridge, 1958).
  * [10] C.H. Llewellyn Smith, Oxford University preprint, to be published in the Proceedings of the 1984 PANIC Conference, Heidelberg.
  * [11] See, for example, the lectures by Gross and by Ellis and the book by Itzykson and Zuber referenced in the preface.
  * [12] The methods employed here are those of light-cone current algebra. More detailed accounts can be found in H. Fritzsch and M Gell-Mann in Broken Scale Invariance and the Light Cone, 1971 Coral Gables Conference on Fundamental Interactions at High Energy, M Dal Cin, G. J. Iverson and A. Perlmutter, eds., (Gordon and Brearh, New York, 1971); see also the lectures by Ellis referenced in the preface.
  * [13] R.L. Jaffe, Nucl. Phys. B229, 205 (1983).
  * [14] The books by Close and Feynman and lectures by Ellis have extensive discussions of the parton phenomenology. For recent experimental summaries see F. Dydak, in Proceedings of the 1983 International Symposium on Lepton and Photon Interactions at High Energies (Newman Laboratories, Ithaca, 1983), D.G. Cassel and D.L. Kreinich, eds., p.634.
  * [15] A. Bodek and A. Simon, University of Rochester preprint UR-906, April 1985, to be published in Z. Phys.
  * [16] Some aspects of the EMC data, notably an enhancement at low$-x$ were anticipated by A. Krzywicki, Phys. Rev. D14, 152 (1976) and R.M. Godbole and D.V.L. Sharma, Phys. Rev. D25, 120 (1982). The majority of papers prior to 1982 treating inelastic electron scattering from nuclei emphasized shadowing ($S^{A}<1$) at low$-x$, an effcct which has yet to be observed at large $Q^{2}$.
  * [17] A. Bodek et al., Phys. Rev. Lett. 50, 1431 (1983), 51, 534 (1983).
  * [18] R. Arnold et al., Phys. Rev. Lett. 52, 727 (1984).
  * [19] A. Argento (BCDMS Collaboration) talk presented at the XX Rencontre de Moriond, Les Ares, 1985.
  * [20] A.M. Cooper et al., Phys. Lett. B 141, 194 (1984). For further references, see the review by A. Bodek in Neutrino 84, University of Rochester Preprint UR-884, June 1984.
  * [21] J.P. Merlo (preliminary) data presented at the XX Rencontre de Moriond, Les Ares, 1985.
  * [22] R.L. Jaffe, Phys. Rev. Lett., 50, 228 (1983).
  * [23] R.L. Jaffe, F.E. Close, R.G. Roberts and G.G. Ross, Phys. Lett. B134, 449 (1984).
  * [24] L.I. Frankfurt and M.I. Strikman, LNPI preprint 886 (1983).
  * [25] Recently, West has argued on rather general grounds at truly asymptotic $Q^{2}$ ($\log Q^{2}/\Lambda^{2}\rightarrow\infty$ where $\Lambda$ is the QCD scale parameter) in QCD, $\delta\epsilon_{\rm Fe}$ should vanish. See G. B. West, Phys. Lett. 54, 2576 (1985). [Such astronomical values of $Q^{2}$ are not obtained experimentally, e.g., $F_{2}(x)=0$ for all $x\neq 0$ in this limit, but West uses the structure of the corrections to this limit to argue rather convincingly that the EMC data should be renormalized downward at low-$x$.]
  * [26] The origins of the convolution model are obscure. It has been well-known to workers in the field for many years. It has been discussed recently by G. B. West, Los Alamos preprint, LA-UR-84-241, C. H. Llewellyn Smith, Phys. Lett. B128, 112 (1983), in Refs. [22] and [23] and in R. L. Jaffe, Comments Nucl. Part. Phys. 13, 39 (1984).
  * [27] G. B. West, Ann. Phys. (NY) 74, 646 (1972); W. B. Atwood and G. B. West, Phys. Rev. D7, 773 (1973); L. I. Frankfurt and M. I. Strickman, it Phys. Lett., 64B, 433 (1976); 65B, 151 (1976); 76B, 333 (1978); Nucl. Phys. B148, 107 (1979); Phys. Lett. 83B, 407 ( 10’/!I), 94B, 216 (1980); A. Bodek and J. L. Ritchie, Phys. Rev. D23, 1070 (1981), D24, 1400 (1981).
  * [28] C.H. Llewellyn Smith, Phys. Lett. 128B, 107 (1983);M. Erickson and A.W. Thomas, Phys. Lett. 128B, 112 (1983); E.L. Berger, F. Coester and R.B. Wiringa, Phys. Rev. D29, 398 (1984).
  * [29] C.E. Carlson and T.J. Havens, Phys. Rev. Lett., 51, 261 (1983), for further references and some discussion see N.N. Nikolaev, Oxford University preprint OT-TP-58/84 (1984).
  * [30] K. Heller and J. Szwed, Acta. Phys. Pol. B16, 157 (1985); H. Faissner and B. R. Kim, Phys. Lett. 130B, 321 (1983).
  * [31] M. Erickson and A. W. Thomas, reference [28].
  * [32] B. L. Birbrair, A. B. Gridnev, M. B. Zhalov, E. M. Levin and V. E. Starodubsky, LNPI preprint 1031 (1985); S. V. Akulinichev, S. A Kulagin and G. M. Vagradov, NBI preprints NBl-85-08 and NBl-85 20 and references therein.
  * [33] For a summary of the equivalent photon approximation see, for example, J.D. Jackson, Classical Electrodynamics, 2nd edition (J. Wiley and Sons, New York, 1975).
  * [34] G. Altarelli and G. Parisi, Nucl. Phys. B126, 298 (1977).
  * [35] M.S. Chen and P. Zerwas, Phys. Rev. D12, 187 (1975).
  * [36] J.B. Kogut and L. Susskind, Phys. Rep. 8C, 77 (1973).
  * [37] Near $x=1$ the radiative process $q\rightarrow qg$ becomes indistinguishable from the elastic process $q\rightarrow q$. To order $\alpha_{c}$ one must include the interference between the gluon vertex correction to Fig. 20 (a) and the zeroth order coupling of the probe to the quarks. As usual with such infrared singularities, the divergence in the vertex correction cancels the divergence in the emission amplitude yielding a finite result.
  * [38] For a discussion see G.B. West, Phys. Rev. Lett. 54, 2576 (1985).
  * [39] M. Gell-Mann and F.E. Low, Phys. Rev. 95, 1300 (1954).
  * [40] For a general analysis of the running coupling and its relation to deep inelastic scattering, see N. Christ, B. Hasslacher and A. H. Mueller, Phys. Rev. D6, 3543 (1973), see also Ref. [36] and references quoted in the preface.
  * [41] It has been shown that the only renormalizable field theories in 4 dimensional space-time that are asymptotically free are non-Abelian theories (see S. Coleman and D.J. Gross, Phys. Rev. Lett. 31, 851 (1973)).
  * [42] K.G. Wilson, Phys. Rev. 179, 1499 (1969). The approach described here is now well-treated in textbooks. See, for example, Cheng and Li and Itzykson and Zuber referenced in the preface.
  * [43] See, for example, R.L. Jaffe and M. Soldate, Phys. Rev. D26, 49 (1982).
  * [44] $O^{\mu_{1}...\mu_{n}}$ actually form an irreducible representation of the Lorentz group labelled $(n/2,n/2)$. It contains representations of the rotation sub-group up to a maximum spin-$n$.
  * [45] S. Weinberg, Phys. Rev. 118, 838 (1960).
  * [46] Operators of odd twist (e.g. $\overline{\psi}\psi$) are not invariant under the chiral transformation $q\rightarrow\exp(i\alpha\gamma_{5})q$ which is a symmetry of QCD with zero mass quarks, so they appear in OPE multiplied by quark masses.
  * [47] R.L. Jaffe and M. Soldate in Proceedings of the 1981 Tallahassee Conference on Perturbative QCD, D. Duke and J. Owens, eds. (AIP, New York, 1981), and [43]; E.V. Shuryak and A. Vainshteyn, Phys. Lett. 105B, 65 (1981); R.K. Ellis, V. Furmanski and R. Petronzio, Nucl. Phys. B207, 1 (1982); B212, 29 (1983).
  * [48] Color indices have been suppressed in Eq. (5.37). $\psi_{a}$ is a three component complex vector in color: $\psi_{ai}$. $D_{\mu}$ is a $3\times 3$ matrix $(D_{\mu})_{ij}=\delta_{ij}\partial_{\mu}-ig\sum_{a=1}^{8}\frac{\lambda^{a}_{ij}}{2}A^{a}_{ij}$ where $A^{a}_{ij}$ is the four vector potential for the gluon field. It is easy to see that only the $+$ component of $O_{n,a}$, i.e., $O_{n,a}^{(\mu^{2})++..+}$ determines $\Theta_{n,a}(\mu^{2})$. It is conventional to use “light-cone gauge”, $A^{+}=0$ so $D^{+}\rightarrow\partial^{+}$.
  * [49] O. Nachtmann, Nucl. Phys. B63, 237 (1973), B78, 455 (1974); R. Barbieri, J. Ellis, M. K. Gaillard and G. G. Ross, Phys. Lett. 64B, 171 (1976), Nucl. Phys. B117, 50 (1976); H. Georgi and H. D. Politzer, Phys. Rev. Lett. 36, 1281 (1976), 37, 68E (1976); A. DeRujula, H. Georgi and H. D. Politzer, Phys. Lett. 64B, 428 (1977), Ann. Phys. (NY) 103, 315 (1977); V. Baluni and E. Eichten, Phys. Rev. Lett. 37, 1181 (1976), Phys. Rev. D14, 3045 (1976).
  * [50] Target mass corrections enter the operator product expansion when the operators are made traceless and symmetric. If the complete tensor in Eq. (5.31) is retained through the analysis, the polynomial in $M^{2}/q^{2}$ which results generates the target mass corrections.
  * [51] These are “current quark” masses. The constituent masses which are larger and appear in non-relativistic quark models are non-perturbative in origin and appear here among higher twist effects.
  * [52] R. L. Jaffe and G. G. Ross, Phys. Lett. 93D, 313 (1980).
  * [53] R. L. Jaffe, Phys. Rev. D11, 1953 (1975).
  * [54] L. Baulieu and C. Kounnas, Nucl. Phys. B155, 429 (1979).
  * [55] F. E. Close, R. G. Roberts and G. G. Ross, Phys. Lett. 120B, 346 (1983).
  * [56] A superficially similar but fundamentally quite different approach was proposed by O. Nachtmann and H.J. Pirner, Z. Phys. C21, 277 (1984).
  * [57] R. L. Jaffe, F. E. Close, R. G. Roberts and G. G. Ross, Phys. Lett. 134B, 449 (1984).
  * [58] F. E. Close, R. L. Jaffe, R. G. Roberts and G. G. Ross, Phys. Rev. D31, 1004 (1985).
  * [59] T. H. Hanssen and R. L. Jaffe, Phys. Rev. D28, 882 (1983).
  * [60] J. Negele, private communication.
  * [61] V. Gribov, in Proceedings $4^{th}$ Winter Seminar on the Theory of the Nucleus and the Physics of High Energies (Academy of Sciences, U.S.S.R., 1969); J. D. Bjorken, in Particle Physics (Irvine Conference 1971), M. Bander, G. L. Shaw and D. Y. Wong, eds., (AIP Conference Proceedngs No. 6, AIP, New York, 1972); S. J. Brodsky in Quarks and Nuclear Forces, Springer Tracts in Modern Physics, No. 100 (Springer-Verlag, Berlin, 1977); V. I. Zakharov and N. N. Nikolaev, Yad. Fiz. 21, 434 (1975), Sov. J. Nucl. Phys. 21, 227 (1975); A.H. Mueller, in Proceedings of the XVII Recontre de Moriond, J. Tran Thanh Van, ed., (Editions Frontières, Dreux, 1982); J.D. Bjorken, (unpublished, private communication from A.H. Mueller).
  * [62] M.S. Goodman, et al., Phys. Rev. Lett. 47, 293 (1981).
  * [63] This description can be formulated in an infinite momentum frame. There the virtual photon does not create pairs but instead scattering from a sea of quark-antiquark pairs. See, for example, the paper of Mueller under [61] above for a discussion of the relation between the two pictures.
  * [64] The following material is derived largely from discussion with A.H. Mueller and unpublished notes by J.D. Bjorken.