# Periodic versus Intermittent Adaptive Cycles in Quasispecies Coevolution

Alexander Seeholzer  Present address: Laboratory of Computational
Neuroscience, EPF Lausanne, 1015 Lausanne, Switzerland.  Arnold-Sommerfeld-
Center für Theoretische Physik and Center for NanoScience Ludwig-Maximilians-
Universität München, Theresienstrasse 37, 80333 München, Germany  Erwin Frey
Arnold-Sommerfeld-Center für Theoretische Physik and Center for NanoScience
Ludwig-Maximilians-Universität München, Theresienstrasse 37, 80333 München,
Germany  Benedikt Obermayer<EMAIL_ADDRESS>Present address:
Max-Delbrück-Center for Molecular Medicine, 13092 Berlin, Germany.  Arnold-
Sommerfeld-Center für Theoretische Physik and Center for NanoScience Ludwig-
Maximilians-Universität München, Theresienstrasse 37, 80333 München, Germany

###### Abstract

We study an abstract model for the coevolution between mutating viruses and
the adaptive immune system. In sequence space, these two populations are
localized around transiently dominant strains. Delocalization or error
thresholds exhibit a novel interdependence because immune response is
conditional on the viral attack. An evolutionary chase is induced by
stochastic fluctuations and can occur via periodic or intermittent cycles.
Using simulations and stochastic analysis, we show how the transition between
these two dynamic regimes depends on mutation rate, immune response, and
population size.

Evolution is commonly pictured as a dynamic process on a fitness landscape in
sequence space. In general, this landscape depends not only on the genotype
but varies dynamically as a function of the environment and coevolving
interaction partners Wright1932 . Prominent biological examples are the
coevolutionary dynamics between the adaptive immune system and virus
populations such as HIV Nowak1990 ; Woo2012  or influenza Bedford2011 , or
between bacteria and their phages Levin2013 . Continuous evolutionary
innovations allow the virus to transiently escape immune suppression,
triggering subsequent adaptations of the immune system. These dynamics can
lead to coevolutionary cycles, which have been generally described in two
different forms Woolhouse2002 ; Bedford2011 : either as an intermittent series
of quasistationary states connected by stochastic jumps, or as periodic and
largely deterministic oscillations. From a modeling perspective, this highly
complex process is determined by three main features Perelson2002 . First,
mutation rates are high and populations are large, which implies large genetic
heterogeneity within the populations Domingo1985 . This has often been
pictured in terms of broad quasispecies distributions around peaks in the
fitness landscape Wang2006 ; Eigen1977 . At the same time, continuous adaption
and coevolutionary arms races are driven by strong ecological interactions
Woolhouse2002 ; Tellier2013 . These modulate effective fitness landscapes
Gavrilets1998 ; Nowak1990 ; Arnoldt2012  and lead to nontrivial nonlinear
population dynamics. Finally, stochastic effects in finite populations become
especially pronounced whenever the first two issues are relevant at the same
time Tsimring1996 ; Rouzine2001 ; Tellier2013 ; Gokhale2014 .

Figure 1: (a) Schematic model for the coevolutionary dynamics of virus (V) and
immune system (IS). The two populations are subject to mutation and selection
(left), but also to ecological interactions (right). (b) Exemplary
trajectories of the relative frequencies of virulent strains $x_{i}$ (top) and
corresponding antibodies $y_{i}$ (bottom). Regular oscillatory dynamics
involving three strains turn into simpler two-strain oscillations at $T_{3\to
2}$ and finally transition into intermittency at $T_{\text{int}}$. Genetic
variability within the populations is calculated from the average pairwise
Hamming distance and indicated in gray dashed lines. (c) Sketch of the
dynamics of the full population distributions in sequence space as described
in the text.

Here, we offer a synthetic perspective on these processes. In our model [see
Fig. 1(a)], we consider a population of $N$ viruses represented by their
genotypes (binary sequences of length $L$ and frequency $x_{i}$) and
replication rates $r_{i}=1$. A small number $n$ of these genotypes
corresponding to particularly virulent strains have a fitness advantage
$\alpha$ over the unit baseline, giving $r_{i}=1+\alpha$ for $i=p,q,\ldots$.
Offspring sequences undergo mutations with per-base rate $\mu_{x}$. In the
absence of immune suppression, and in the stationary state, the viral
population localizes as so-called quasispecies around any of the fittest
genotypes, provided the mutation rate is smaller than Eigen’s error threshold
$\mu_{\text{c}}\approx\ln(\alpha+1)/L$ Eigen1977 . This simple picture is
considerably complicated by the host’s adaptive immune system, which produces
antibodies that recognize and neutralize viruses with matching epitopes
Alberts2002 . Antibody production is specifically increased and variability in
the binding affinity is introduced when viruses with matching genotype are
encountered Rajewsky1996 , in a process that can be modeled in terms of
mutation and selection. Similar concepts can be used for bacterial immune
systems, where spacer sequences in the host genome complementary to genetic
elements of a phage take antibodylike functions Levin2013 . Hence, for the
immune system we introduce a second population of $N$ binary sequences with
frequencies $y_{i}$, mutation rate $\mu_{y}$, unit replication rate for
unstimulated production, and stimulated antibody production in the presence of
perfectly matching viruses Kamp2002 ; Wang2006 . Ecological interactions then
introduce frequency-dependent fitness terms $\propto x_{i}y_{i}$ for such
matched virus and antibody pairs. Including these terms leads to a reduction
of the viral load and stimulation of antibody production [Fig. 1(a), right].
In the deterministic limit ($N\to\infty$), our model is described by
supplement

$\begin{split}\dot{x}_{i}&=\textstyle{\sum_{j}}m^{x}_{ij}r_{j}x_{j}-\alpha
x_{i}y_{i}-x_{i}\phi_{x},\\\
\dot{y}_{i}&=\textstyle{\sum_{j}}m^{y}_{ij}(1+\gamma
x_{j})y_{j}-y_{i}\phi_{y},\end{split}$ (1)

where $i$ and $j$ run over all $2^{L}$ sequences. The fitness advantage
$\alpha$ of virulent strains can be suppressed to background levels by a
perfectly adapted immune system, which undergoes stimulated production at rate
$\gamma$ when encountering matching viral epitopes. Further,
$m^{x}_{ij}=\mu_{x}^{d_{ij}}(1-\mu_{x})^{L-d_{ij}}$ is the probability of
having $d_{ij}$ simultaneous mutations, where $d_{ij}$ is the Hamming distance
between $x_{i}$ and $x_{j}$. The dilution terms $\phi_{x/y}$ are obtained from
the conditions $\sum_{i}\dot{x}_{i}=0$ and $\sum_{i}\dot{y}_{i}=0$,
respectively, and keep the sizes of the two populations fluctuating around
constant values. This constraint applies to the stationary phase of the
adaptive race, while we ignore some of the effects of a changing viral load
Nowak1990 ; Woo2012 ; Wang2006 ; Bull2007  and also neglect immune system
memory Bianconi2010  and unspecific recognition Alberts2002 .

To facilitate a systematic study of the effects of demographic noise by means
of simulations and theoretical analysis, our starting point is the underlying
stochastic master equation supplement  which has rarely been used in this
context. Its deterministic limit leads to Eq. (1) and connects to established
quasispecies theory Eigen1977 ; Kamp2002 . Exemplary simulation results
obtained with the Gillespie algorithm Gillespie1992  are shown in Fig. 1(b),
with parameters in the coexistence regime discussed below. We readily identify
characteristics of the intermittent coevolutionary dynamics. First, a
particularly virulent strain with its associated quasispecies “cloud” of
mutants triggers a specific immune response (a), leading to a corresponding
localization in the antibody sequence space (b). This gives alternative viral
strains that are not under immune attack a fitness advantage, and after a
brief “search” period during which the viral population becomes delocalized,
this new fitness peak is colonized in a “growth” phase (c), awaiting the
adaptive immune response (d). The delocalization and relocalization dynamics
of each population in sequence space are clearly visible as transient
increases in their respective mean pairwise Hamming distances [Fig. 1(b)].
Intriguingly, this sequence of events can occur both in the form of regular
oscillations as well as by means of stochastically intermittent cycles
Woolhouse2002 . The former occurs when the large genetic diversity within the
population extends across the valleys between different fitness peaks and
signifies periodic shifts in the extent to which these peaks are populated
Gavrilets1998 . The latter case indicates that adaptation proceeds
stochastically via the random discovery of previously unpopulated fitness
peaks by relatively tightly localized populations.

## Steady-state regimes: coexistence for mutation rates below interdependent
error thresholds.

We use a reduced deterministic version of the model to determine stationary
states and the associated error thresholds. We restrict the analysis to the
populations of the $n$ virulent strains $x_{p,q,\ldots}$ and their respective
antibodies $y_{p,q,\ldots}$, and lump all mutant sequences together in the so-
called error tail Schuster1988 . The high-dimensional system (1) is then
reduced to supplement

$\begin{split}\dot{x}_{p}&=\left[Q_{x}\left(1+\alpha\right)-\alpha
y_{p}-\bar{\phi}_{x}\right]x_{p},\\\ \dot{y}_{p}&=\left[Q_{y}\left(1+\gamma
x_{p}\right)-\bar{\phi}_{y}\right]y_{p},\end{split}$ (2)

where $p$ runs over the $n$ strains which are coupled by the corresponding
dilution terms $\bar{\phi}_{x/y}$. $Q_{x/y}=(1-\mu_{x/y})^{L}$ are the quality
factors. A straightforward stability analysis of fixed points in this system
with respect to $\mu_{x/y}$ as bifurcation parameters yields the phase
diagrams of Fig. 2.

Figure 2: (a) Regimes of coevolution. High mutation rates $\mu_{x}$ of the
virus lead to population delocalization, while for lower mutation rates a
regime of coexistence emerges. Intermediate values lead to a degenerate
localization regime for the virus (see text). (c) Steady-state values for
relative frequencies $x_{p}=x_{q}$ and $y_{p}=y_{q}$ as a function of
$\mu_{y}$ with $\mu_{x}=\mu_{y}$ (above) or $\mu_{x}=0.05$ (below). Solid
lines are solutions of Eq. (2) and dots are simulation results. Panels (a) and
(c) are for $\gamma=\alpha$, while (b) and (d) show analogous result for
$\gamma=10\alpha$, where $L=8$, $n=2$, and $\alpha=10$.

As expected, we recover the classical result that the viral population
localizes around a fitness peak only if
$Q_{x}>Q_{c}\equiv\left(\alpha+1\right)^{-1}$, with increasing genetic
variability (i.e., the width of the population distribution) for larger
mutation rate $\mu_{x}$. However, antibodies are localized only (1) if their
mutation rate $\mu_{y}$ is small enough, (2) if their production rate $\gamma$
is high enough and (3) if the virus attack is specific enough (i.e., tightly
localized). These interdependent requirements are an inevitable consequence of
ecological interactions, and they translate into the condition
$Q_{y}=\left\\{\left(\gamma/\alpha
n\right)\left[Q_{x}(\alpha+1)-1\right]+1\right\\}^{-1}$ as the analytical
limit for the coexistence regime (blue dashed lines in Fig. 2). Only in this
regime do we find the intriguing oscillatory dynamics shown in Fig. 1 that
will be discussed below. Finally, in a somewhat model-specific “degenerate”
regime bounded by
$Q_{y}=\left\\{\left(\gamma/\alpha\right)\left[Q_{x}(\alpha+1)-1\right]+1\right\\}^{-1}$,
the virus population can stably localize about several fitness peaks
simultaneously such that none of these quasispecies is sufficiently tight to
trigger a specific response of the immune system, which thus remains
delocalized. The fixed points of the approximate system (2) coincide closely
with the mean steady-state concentrations obtained by stochastic simulation of
the full system (1) [see Figs. 2(c) and 2(d)]. Interestingly, for
$\alpha=\gamma$ and symmetric mutation rates $\mu_{x,y}\equiv\mu$, the
critical condition of coexistence can be approximated by
$\mu\approx\left(1/2L\right)\ln\left(\alpha/2\right)$ for large $\alpha$ and
$L$, which generalizes a comparable result for mutualistic frequency-dependent
fitness Obermayer2009  to the case of antagonistic interactions. This
correspondence also suggests that the error thresholds derived here should be
largely unchanged if recognition between the two population tolerates some
mismatches Obermayer2010 .

## Noise-driven oscillations in the coexistence regime.

Performing a linear stability analysis in the coexistence regime reveals that
the oscillations seen in the simulations are caused by $n-1$ pairs of purely
imaginary eigenvalues. Numerical solutions of the deterministic Eqs. (2) show
complex but regular oscillations involving all $n$ strains with slow amplitude
variations controlled by higher-order nonlinearities (see Fig. S1 in the
Supplemental Material supplement ). Results from stochastic simulations,
however, suggest that such complex patterns quickly give rise to simpler
oscillations involving only two strains, which at a later time transition into
intermittency (cf. Fig. 1). Investigating the case $n>2$ by simulations below,
we restrict further analysis to $n=2$. Also, here we only display more compact
analytical results for the case $\gamma=\alpha$ (see the Supplemental Material
supplement  for general results). Our analysis exploits that in the
coexistence regime mutation rates $\mu_{x/y}\lesssim\ln\alpha/L$ are small
compared to the error thresholds and can be used as expansion parameters. To
obtain the nonlinearities that control oscillation amplitudes, we expand Eq.
(2) to first order and transform to polar normal form on the two-dimensional
stable manifold supplement :

$\displaystyle\dot{u}$
$\displaystyle=-\frac{4}{5}L\left[\mu_{x}\left(\alpha+1\right)+\mu_{y}\right]u^{2},$
(3a) $\displaystyle\dot{\varphi}$
$\displaystyle=\frac{\alpha}{2}-\frac{L}{2}(\alpha+1)(\mu_{x}+\mu_{y})+\mathcal{O}(u),$
(3b)

where $u$ is a squared radial coordinate indicating deviations from the
coexistence fixed point and $\varphi$ measures the phase of the oscillations.
Equation (3a) exhibits a weak geometric decay of the oscillation amplitude
$\mathcal{O}(u^{2}\mu_{x})\ll 1$, which makes the fixed point only marginally
stable and thus vulnerable to stochastic fluctuations Reichenbach2006 ;
Cremer2009 ; Bladon2010 . Notably, the oscillation frequency of Eq. (3b)
depends mainly on the fitness advantage $\alpha$, and is only weakly slowed
down by mutations. In this deterministic regime, the quasispecies distribution
in sequence space is broad enough that the time required to shift to a new
fitness peak is dominated by the _growth_ of the subpopulation already on the
new peak (with a rate $\alpha$) rather than the _search_ for this new peak in
the first place (via mutations). We note that this effect is even stronger if
the two fitness peaks are close in sequence space, i.e., if direct mutations
between them are not ignored as in Eq. (2). In contrast, when the coexistence
regime displays intermittent dynamics, because the relevant sequence space is
not already inhabited by the virus population, the dynamics are inherently
stochastic and mutation rates can be _too small_ for the virus to explore
enough sequence space to escape immune suppression in time. This would
correspond to an adaptation threshold as found in a previous study Kamp2002 .
However, as shown more formally below, this situation is incompatible with the
presence of deterministic dynamics, which is an assumption of standard
quasispecies theory. Instead, population genetics models should be used
Tellier2013 ; Arnoldt2012 .

## Noise determines if dynamics are periodic or intermittent.

We can characterize how stochastic noise controls the transition between
periodic and intermittent adaptive dynamics by means of stochastic averaging.
This technique enables a systematic derivation of effective one-dimensional
Fokker-Planck equations in relevant subspaces of more complex high-dimensional
nonlinear dynamics such as those arising in evolutionary game theory
Dobrinevski2012 . It is based on the time scale separation between slow radial
and fast azimuthal dynamics in Eq. (3): $\varphi$ evolves on fast time scales
($\dot{\varphi}\propto\alpha$), while $u$ changes much more slowly
($\dot{u}\propto\mu Lu^{2}$). Using this observation, we can derive effective
coefficients governing the evolution of the probability distribution $P(u,t)$
of the radial variable by averaging the angular dynamics over one oscillation
period supplement . To leading order, we get

$\partial_{t}P=-\partial_{u}\left[\left(-a_{1}u^{2}+\frac{a_{2}}{N}\right)P\right]+\frac{1}{N}\partial_{u}^{2}\left(a_{2}u\,P\right),$
(4)

with $a_{1}=\tfrac{4}{5}L[\mu_{x}\left(\alpha+1\right)+\mu_{y}]$ and
$a_{2}=\frac{1}{16}\big{\\{}4+3\alpha-\mu_{x}L\left[\left(4/\alpha\right)+11+7\alpha\right]-\mu_{y}L\left[\left(4/\alpha\right)+7+3\alpha\right]\big{\\}}$.
Note that in the deterministic limit $N\to\infty$ we recover Eq. (3a). For a
finite population, we now find the deterministic decay $\propto a_{1}u^{2}$
towards the coexistence fixed point in competition with a stochastic outward
drift $\propto\left(a_{2}/N\right)$, which destabilizes the fixed point and
leads to a finite oscillation amplitude $\langle
u\rangle=\sqrt{\left(2/\pi\right)\left(a_{2}/a_{1}N\right)}$ . As mutation
rates get small, expected oscillation amplitudes grow as
$[(\alpha+1)\mu_{x}+4\mu_{y}]^{-1/2}$, eventually hitting the borders of the
concentration simplex. This indicates the transition from regular oscillations
to intermittent behavior: during large-amplitude oscillations the fittest
virus genotypes are temporarily lost from the population and are only much
later recovered through spontaneous mutations.

A more detailed understanding of this transition is obtained by estimating the
lifetime of the regular oscillations. To this end, we use the bounds on the
radial variable $u_{\text{max}}=\tfrac{1}{8}-\mathcal{O}(\mu L)$, where the
populations are fully localized about only one peak. The chances of observing
a transition to intermittent behavior are estimated from the mean first
passage time (MFPT) $T_{\text{int}}$ from $u=0$ to $u=u_{\text{max}}$ under
Eq. (4). Using standard methods Kampen1992 , we find the result supplement

$T_{\text{int}}=N\frac{u_{\text{max}}}{a_{2}}\tilde{F}\left(\frac{Na_{1}u^{2}_{\text{max}}}{2a_{2}}\right),$
(5)

where ${\tilde{F}}(x)$ is the generalized hypergeometric function
${{}_{2}F_{2}}(\nicefrac{{1}}{{2}},1;\nicefrac{{3}}{{2}},\nicefrac{{3}}{{2}};x)$.
Equation (5) can be brought into scaling form by defining
$N^{*}=\left(2a_{2}/a_{1}u^{2}_{\text{max}}\right)$ and
$T^{*}=N^{*}\left(u_{\text{max}}/a_{2}\right)$. To compare this result to
simulations, we plot the rescaled MFPT
$\left(T_{\text{int}}/T^{*}\right)=\left(N/N^{*}\right){\tilde{F}}(N/N^{*})$
(see Fig. 3). The nearly perfect data collapse for different parameter choices
validates our analytical approach.

While $N^{*}$ measures the population size at the crossover from periodic to
intermittent dynamics, $T^{*}$ denotes the corresponding typical duration of
the transition. For large populations ($N>N^{*}$), we find $T_{\text{int}}\sim
N^{-1/2}e^{N/N^{*}}$; this almost exponential growth of the MFPT indicates
that the dynamics are effectively deterministic and intermittent behavior
extremely unlikely. For $N<N^{*}$, we find $T_{\text{int}}\sim N$ and the
dynamics thus easily transition into intermittency. This distinction based on
the scaling of $T_{\text{int}}$ with $N$ has recently been suggested in the
context of game theory Cremer2009 . In our case, however, finite mutation
rates prevent permanent extinction of subpopulations and stabilize regular
oscillatory behavior even in small populations, because the deterministic
decay in Eq. (3a) is strengthened and the critical population size
$N^{*}\propto(\mu_{x}L)^{-1}$ is reduced. Thus, even for small populations,
mutations can act as a driving force for the stabilization of regular
oscillations, which _a posteriori_ justifies assumptions underlying
quasispecies theory and generalizes previous observations Gokhale2014 . In
contrast, from results for general $\gamma$ (see Fig. S2 in the Supplemental
Material supplement ), we find that a strong immune response (i.e.,
$\gamma>\alpha$) promotes early transitions into intermittency [cf. Fig.
3(b)], since both $N^{*}$ and $T^{*}$ increase with $\gamma$. However, these
parameters are insensitive against the precise value of $\mu_{y}$ [cf. Fig.
3(a)]; this suggests that effective immune suppression is achieved via a
strong stimulated response rather than high adaptive flexibility. Indeed,
extreme antibody secretion rates have been reported in the literature Hibi1986
. Finally, we support our choice of limiting the analysis to $n=2$ strains by
simulating a system with $n=3$ strains, measuring the time $T_{3\to 2}$ until
one strain is lost as well as the subsequent $T_{\text{int}}$ until the
remaining two strains transition to intermittency. As shown in Fig. 3(c), the
state with all three strains present is short lived compared to the two-strain
oscillations, especially in the relevant deterministic regime of larger
population size. Hence, apart from numerical prefactors the general trend
captured in Eq. (5) also describes systems with larger $n$.

Figure 3: Mean time until transition from regular oscillations to
intermittency. Dashed lines are from the analytical result (5). (a) Rescaled
simulation data for $L=8$, $\alpha=\gamma=10$ and different choices of
$\mu_{x/y}$ collapse onto a universal curve (unscaled data shown in the
inset). (b) Transition times decrease for increasing $\gamma$ [parameters
otherwise as in (a)]. (c) For $n=3$ virulent strains, the transition time
$T_{3\to 2}$ until one strain is lost is much shorter than $T_{\text{int}}$,
especially for large populations.

## Conclusions.

We have analyzed a model for the coevolutionary dynamics of virus and immune
system, combining simulations with nonlinear deterministic and stochastic
analysis. Starting from the established quasispecies treatment of this
problem, we explicitly introduced interactions between the populations. These
lead to interdependent error thresholds, because a focused immune defense
against a specific viral strain is impaired for large genetic variability in
the virus population. Further, we performed a rigorous analysis of stochastic
effects in the coexistence regime: regular yet noise-induced oscillatory
behavior for large populations, large mutation rates, and weak immune response
turn into stochastic intermittent cycles for smaller populations, smaller
mutation rates and strong immune response. Our simulations indicate that the
reverse transition from intermittency towards regular oscillations is a rare
event occurring on time scales well beyond $T_{\text{int}}$. It cannot easily
be analyzed within our reduced two-dimensional model as it will depend on the
entire population structure. Finally, we note that our abstract model based on
quasispecies theory focuses on the dynamics of genetic variability within
populations of constant size. This assumption is of course violated for some
biological scenarios, where immune response modulates the viral load Nowak1990
; Wang2006 ; Bull2007  and may well lead to extinction of the virus Tejero2010
; Woo2012 . We expect that more detailed models including these and other
effects relevant in biological situations Bianconi2010 ; Alberts2002  will
also be amenable to theoretical analysis based on the stochastic averaging
techniques used here.

###### Acknowledgements.

We acknowledge helpful comments by anonymous reviewers and financial support
by the Deutsche Forschungsgemeinschaft in the framework of the SFB/TR12 –
Symmetries and Universality in Mesoscopic Systems.

## References

  * (1) S. Wright, in Proceedings of The Sixth International Congress of Genetics, Ithaca, New York, 1932 (George Banta Publishing, Menasha, WI, 1932), 356.
  * (2) M. A. Nowak, R. May, and R. Anderson, AIDS 4, 1095 (1990).
  * (3) H.-J. Woo and J. Reifman, Proc. Natl. Acad. Sci. U.S.A. 109, 12980 (2012).
  * (4) T. Bedford, S. Cobey, and M. Pascual, BMC Evol. Biol. 11, 220 (2011)
  * (5) B. R. Levin, S. Moineau, M. Bushman, and R. Barrangou, PLoS Genet. 9, e1003312 (2013)
  * (6) M. E. J. Woolhouse, J. P. Webster, E. Domingo, B. Charlesworth, and B. R. Levin, Nat. Genet. 32, 569 (2002).
  * (7) A. S. Perelson, Nat. Rev. Immunol. 2, 28 (2002).
  * (8) E. Domingo et al., Gene 40, 1 (1985).
  * (9) G. Wang and M. W. Deem, Phys. Rev. Lett. 97, 188106 (2006).
  * (10) P. Schuster and M. Eigen, Naturwissenschaften 64, 541 (1977).
  * (11) A. Tellier, S. Moreno-Gamez, and W. Stephan, Evolution 68, 2211 (2014).
  * (12) S. Gavrilets and A. Hastings, J. Theor. Biol. 191, 415 (1998).
  * (13) H. Arnoldt, M. Timme, and S. Grosskinsky, J. R. Soc. Interface 9, 3387 (2012).
  * (14) C. S. Gokhale, A. Papkou, A. Traulsen, and H. Schulenburg, BMC Evol. Biol. 13, 254 (2014).
  * (15) L. S. Tsimring, H. Levine, and D. A. Kessler, Phys. Rev. Lett. 76, 4440 (1996).
  * (16) I. Rouzine, J. Coffin, and A. Rodrigo, Microbiol. Mol. Biol. R. 65, 151 (2001).
  * (17) B. Alberts et al., Molecular Biology of the Cell (Garland Science, New York, 2002).
  * (18) K. Rajewsky, Nature (London) 381, 751 (1996).
  * (19) C. Kamp and S. Bornholdt, Phys. Rev. Lett 88, 068104 (2002).
  * (20) See Supplemental Material for details, which includes Refs. Eigen1978 ; Jones1981 ; Nowak1989 ; Andrade1993 ; Hinch1991 ; Wiggins2003 ; Gardiner2009 .
  * (21) M. Eigen and P. Schuster, Naturwissenschaften 65, 7 (1978).
  * (22) B. Jones and H. Leung, Bull. Math. Biol. 43, 665 (1981).
  * (23) M. A. Nowak and P. Schuster, J. Theor. Biol. 137, 375 (1989).
  * (24) M. Andrade, J. C. Nuño, F. Morán, F. Montero, G. J. Mpitsos, Physica (Amsterdam) 63D, 21 (1993).
  * (25) E. J. Hinch, Perturbation Methods (Cambridge University Press, Cambridge, England, 1991)
  * (26) S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos (Springer, New York, 2003).
  * (27) C. Gardiner, Handbook of Stochastic Methods (Springer, New York, 2009).
  * (28) J. J. Bull, R. Sanjuán, and C. O. Wilke, J. Virol. 81, 2930 (2007).
  * (29) G. Bianconi, D. Fichera, S. Franz, and L. Peliti, J. Stat. Mech. 8, 08022 (2011).
  * (30) D. Gillespie, Physica (Amsterdam) 188A, 404 (1992).
  * (31) P. Schuster and J. Swetina, Bull. Math. Biol. 50, 635 (1988).
  * (32) B. Obermayer and E. Frey, Europhys. Lett. 88, 48006 (2009).
  * (33) B. Obermayer and E. Frey, J. Theor. Biol. 267, 653 (2010).
  * (34) T. Reichenbach, M. Mobilia, and E. Frey, Phys. Rev. E 74, 051907 (2006).
  * (35) J. Cremer, T. Reichenbach, and E. Frey, New J. Phys. 11, 093029 (2009).
  * (36) A. J. Bladon, T. Galla, and A. J. McKane, Phys. Rev. E 81, 066122 (2010).
  * (37) A. Dobrinevski and E. Frey, Phys. Rev. E 85, 051903 (2012).
  * (38) N. G. van Kampen, Stochastic Processes in Physics and Chemistry (Elsevier, New York, 1992).
  * (39) T. Hibi and H. M. Dosch, Eur. J. Immunol. 16, 139 (1986).
  * (40) H. Tejero, A. Marín, and F. Montero, J. Theor. Biol. 262, 733 (2010).

See pages ,1,,2,,3,,4,,5,,6,,7,,8,,9,,10,,11, of supplement_proof.pdf


# Do Machines and Humans Focus on Similar Code? Exploring Explainability of
Large Language Models in Code Summarization

Jiliang Li  Vanderbilt UniversityNashvilleTennesseeUSA
<EMAIL_ADDRESS>,  Yifan Zhang  Vanderbilt
UniversityNashvilleTennesseeUSA<EMAIL_ADDRESS>,  Zachary Karas
Vanderbilt UniversityNashvilleTennesseeUSA<EMAIL_ADDRESS>,  Collin
McMillan  University of Notre DameNotre DameIndianaUSA<EMAIL_ADDRESS>,  Kevin
Leach  Vanderbilt UniversityNashvilleTennesseeUSA<EMAIL_ADDRESS>and  Yu Huang  Vanderbilt UniversityNashvilleTennesseeUSA
<EMAIL_ADDRESS>
(2024)

###### Abstract.

Recent language models have demonstrated proficiency in summarizing source
code. However, as in many other domains of machine learning, language models
of code lack sufficient explainability — informally, we lack a formulaic or
intuitive understanding of what and how models learn from code. Explainability
of language models can be partially provided if, as the models learn to
produce higher-quality code summaries, they also align in deeming the same
code parts important as those identified by human programmers. In this paper,
we report negative results from our investigation of explainability of
language models in code summarization through the lens of human comprehension.
We measure human focus on code using eye-tracking metrics such as fixation
counts and duration in code summarization tasks. To approximate language model
focus, we employ a state-of-the-art model-agnostic, black-box, perturbation-
based approach, SHAP (SHapley Additive exPlanations), to identify which code
tokens influence that generation of summaries. Using these settings, we find
no statistically significant relationship between language models’ focus and
human programmers’ attention. Furthermore, alignment between model and human
foci in this setting does not seem to dictate the quality of the LLM-generated
summaries. Our study highlights an inability to align human focus with SHAP-
based model focus measures. This result calls for future investigation of
multiple open questions for explainable language models for code summarization
and software engineering tasks in general, including the training mechanisms
of language models for code, whether there is an alignment between human and
model attention on code, whether human attention can improve the development
of language models, and what other model focus measures are appropriate for
improving explainability.

Neural Code Summarization, Language Models, Explainable AI, SHAP, Human
Attention, Eye-Tracking

††conference: 46th International Conference on Program Comprehension; April
2024; Lisbon, Portugal††journalyear: 2024††copyright: acmlicensed††conference:
32nd IEEE/ACM International Conference on Program Comprehension; April 15–16,
2024; Lisbon, Portugal††booktitle: 32nd IEEE/ACM International Conference on
Program Comprehension (ICPC ’24), April 15–16, 2024, Lisbon, Portugal††doi:
10.1145/3643916.3644434††isbn: 979-8-4007-0586-1/24/04††ccs: Computing
methodologies Artificial intelligence

##  1\. Introduction

Recent language models for code have shown promising performance on several
code-related tasks (xu2022systematic, ). Among these tasks is neural code
summarization, where a language model generates a short natural language
summary describing a given code snippet. This is often an indicative task
demonstrating a model’s ability to comprehend code. Currently, the majority of
assessments for how well a language model understands code directly measures
the quality of code summaries generated by the models, and compares them with
human-written summaries (xu2022systematic, ). Comparatively little is known
about why and how the language models reason about code to generate such
summaries. Similar to many other downstream domains of machine learning in
software engineering, understanding and explaining how and why language models
for code work (or fail) is critical to improving model architecture, reducing
bias, and preventing undesirable model behavior.

Human programmers typically achieve a strong understanding of code. Thus,
proficient language models might be explained if they focus on the same parts
of code that humans would (paltenghi2021thinking, ). Eye-tracking studies have
been conducted to analyze programmers’ visual patterns while reading code
(abid2019using, ; rodeghero2015empirical, ). Specifically, the duration and
frequency of a programmer’s eye gaze on a part of code in a spatially-stable
manner, referred to as fixation duration and fixation count respectively, are
indicative of cognitive load (sharafi2015eye, ). Thus, these measures of eye-
tracking can indicate the parts of code on which human programmers focus. In
contrast, there is a lack of consensus on how to measure a language model’s
reasoning about code (see Section 2.2). Most existing works extract the self-
attention layers in language models for code to measure the model attention
(paltenghi2021thinking, ; paltenghi2022extracting, ; huber2023look, ). Such
methods require direct access to the internal layers of a language model,
limiting the possibility to investigate interpretability of many state-of-the-
art proprietary models (e.g., ChatGPT).

In this paper, to investigate how proprietary language models reason about
code, we employ a state-of-the-art perturbation-based method, SHAP
(lundberg2017unified, ) (SHapley Additive exPlanations), that treats each
language model as a black-box function. With SHAP, we analyze the feature
attribution (i.e., which parts of code are deemed important by the model) in
six different state-of-the-art language models for code. We use a set of Java
methods to task both the language models and human programmers with writing
code summaries. The feature attribution in the language models, measured by
SHAP, is then compared with human developers’ focus, collected from eye-
tracking. We hypothesize that sufficiently large models may learn to focus on
parts of code similarly to humans. If validated, language model behavior can
thus be described in terms of human behavior, ultimately helping to explain
and improve language models. However, we find that explainability cannot be
provided through this lens and find no statistically significant evidence
suggesting the hypothesized alignment. Furthermore, we did not find that
language models’ focus exhibits a statistically significant correlation with
human focus in general. For future research that aims to explore the
explainability of language models for code summarization, especially for those
leveraging human attention, our findings might suggest the following: (1)
though widely used in AI, SHAP may not be an optimal method to investigate
where language models focus during code summarization, or alternatively, (2) a
misalignment between language models and human developers in reasoning about
code may provide insights for improving AI models for code summarization.

##  2\. Background and Related Work

###  2.1. Neural Models for Code Summarization

Advancements in deep learning have enabled machine learning models to generate
summaries for source code. Among the state-of-the-art models, NeuralCodeSum
(NCS) first introduced the use of Transformers in neural code summarization
(ahmad2020transformer, ). With the rise of large language models (LLMs),
ServiceNow and HuggingFace released a 15.5B parameter LLM for code, StarCoder
(li2023starcoder, ), and Meta released a 7B parameter LLM, Code LLama
(roziere2023code, ), both of which can serve to summarize code. Although not
inherently an LLM for code, GPT3.5 (openai2023chatgpt, ) and GPT4
(openai2023gpt, ) are also capable of code summarization. In this paper, we
investigate how all the aforementioned models reason about code when tasked to
generate code summaries.

###  2.2. Interpretability of Language Models

Existing works on interpretable language models generally seek to investigate
the relative importance of each input token for model performance
(vaswani2017attention, ; hooker2018evaluating, ; molnar2020interpretable, ).
Such works can be commonly categorized into two types: white-box vs. black-
box. White-box approaches require access to a language model’s internal layers
(shrikumar2017learning, ; sundararajan2017axiomatic, ), often directly
investigating the self-attention scores in Transformer-based models
(galassi2020attention, ; zhang2022does, ; zeng2022extensive, ). However,
Transformer-based models’ inherent complexity has led to a lack of consensus
on how to aggregate attention weights (wang2021wheacha, ; zhang2022diet, ;
zhang2022does, ). For the general research community, white-box approaches
preclude proprietary models (e.g., ChatGPT).

In contrast, state-of-the-art black-box approaches like SHAP
(lundberg2017unified, ) (SHapley Additive exPlanations) apply game-theoretic
principles to assess the impact of input variations on a model’s output. SHAP
evaluates the effects of different combinations of input features — such as
tokens in a text sequence — by observing how their presence or absence
(simulated by token masking) alters the model’s prediction from an expected
result. This process helps to ascertain the relative contribution of each
feature to the output, allowing for an analysis of the model without requiring
access to its internal architecture (liu2023reliability, ; kou2023model, ). In
this paper, to investigate proprietary models, we employ SHAP to measure where
language models focus on code.

###  2.3. Comparing Human vs. Machine Attention

Previous papers have examined the alignment between human and model attention
in code comprehension tasks. Paltenghi et al. (paltenghi2022extracting, )
found that CodeGen’s (nijkamp2022codegen, ) self-attention layers attend to
similar parts of code compared to human programmers’ visual fixations when
answering comprehension questions about code. Similarly, Huber et al.
(huber2023look, ) discovered overlaps in attention patterns between neural
models and humans when repairing buggy programs. Notably, Paltenghi and
Prasdel (paltenghi2021thinking, ) compared language models’ self-attention
weights and humans’ visual attention during code summarization. They found
that model attention, measured by self-attention weights, does not align well
with human attention. However, this work is limited by investigating only
small CNN and transformer models. Most importantly, all aforementioned studies
used white-box approaches towards interpretability of open-source models,
limiting applicability to state-of-the-art proprietary models.

Recently, Kou et al. (kou2023model, ) utilized both white-box and black-box
perturbation-based approaches to measure LLMs’ focus in code generation tasks,
and discovered a consistent misalignment with humans’ attention. In general,
these works have demonstrated that whether human and machine attention align
depends heavily on the methods employed to approximate machine focus, as well
as the specific code comprehension task examined. In this paper, we build upon
former works by examining whether human attention correlates with feature
attribution in language models, measured by a black-box perturbation-based
approach, in code summarization.

##  3\. Experimental Design

###  3.1. Measuring Human Visual Focus

We used eye-tracking data measuring human attention from a controlled human
study with 27 programmers. The study obtained IRB approval, and asked
participants to read Java methods and write accompanying summaries
(bansal2023modeling, ). Each participant summarized 24–25 Java methods from
the FunCom dataset (leclair2019recommendations, ), yielding 671 trials of eye-
tracking data in total. Considering data quality, two authors with five and
eight years of Java experience cooperatively removed participant data
associated with five summaries that did not demonstrate an understanding of
the Java code.

In this work, we sought to measure where humans and language models focus on
code as they summarize it. We first used the srcML parser to convert each Java
method into its corresponding Abstract Syntax Tree (AST) representation
(collard2013srcml, ). The AST provides structural context for each token
literal (i.e., ‘Hello World’ $\xrightarrow{}$ String Literal). With the gaze
coordinates collected from the eye-tracker (Tobii_2023, ), we measured humans’
focus on each AST token. Typically, researchers use fixations to quantify
human visual focus (sharafi2015eye, ). A fixation is defined as a spatially
stable eye-movement lasting 100–300ms. Most cognitive processing occurs during
fixations (sharafi2015eye, ), so researchers consider their frequency and
duration in making inferences about human cognition. In our analyses, we
computed the average count and duration of programmers’ fixations on each AST
token. Consequently, for each Java method, we obtained two visual focus
vectors with lengths equal to the number of AST tokens, respectively, which
represent fixation counts and durations on each token111Our analyses do not
include brackets or semi-colons, or other such syntactic elements..

###  3.2. Measuring Model Focus

As mentioned in Section 2.2, we choose SHAP’s official, default implementation
of the TeacherForcing method to measure feature attribution in language
models, treating each as a black-box function. For each language model, we
pass in each of the 68 Java methods (also read by human programmers) as input,
along with necessary prompting for the model to output summaries of source
code. For each Java method passed into each language model, we let $i$ denote
an input token (in code) and $o$ denote an output token (in summary). For each
($i$, $o$) pair, SHAP produces an importance score, denoted $v_{(i,o)}$,
signifying how much $i$’s presence or absence alters the presence of $o$.
Then, the importance score of each input token222We use the absolute value by
choice, without which experiments show similar results., $v_{i}$, is
calculated such that $v_{i}=\sum_{o}|v_{(i,o)}|$. Note that now $v_{i}$ is
associated with a language model token, and each AST token may consist of
several language model tokens. Thus, for each AST token, we calculate its
score $\frac{1}{n}\sum_{j=1}^{n}v_{j}$, where $v_{1},\cdots,v_{n}$ are scores
of language model tokens constituting the AST token. Consequently, for each
language model on each Java method, we obtain a focus vector (with a length
equal to the number of AST tokens) representing how influential each AST token
is to the model.

In total, we investigated the model focus of six different models: GPT4, GPT-
few-shot, GPT3.5, StarCoder, Code Llama, and NCS. Here, GPT-few-shot is a
GPT3.5 model, but in an attempt for the model to produce code summaries more
similar to those of humans, we used few-shot prompting to instruct the model
to provide summaries similar to two randomly selected human-written summaries.
The other five state-of-the-art LLMs are introduced in Section 2.1 and
implemented with their default parameters.

###  3.3. Comparing Human and Model Foci

For brevity, we refer to two human visual focus measurements (i.e., fixation
duration and count) and six language models as eight ”focus sources.” For each
source, we obtained 68 focus vectors, each corresponding to a Java method.
These vectors were normalized to sum to $1$, and reflect how important each
AST token is for the human/model. We answer these research questions:

  * $\bullet$

RQ1: Is there a general correlation between human and machine focus patterns
for code summarization?

  * $\bullet$

RQ2: Do the code summaries increase in quality when machine focus becomes more
aligned with that of humans?

####  3.3.1. RQ1

We assess the correlation between human and machine foci across the 68 Java
methods. Specifically, for each pair of focus sources, we iterate through each
Java method and calculate the Spearman’s rank coefficient ($\rho$)
(spearman1961proof, ) between the two sources’ vectors for that method. Then,
for each pair of focus sources, we report: (1) The mean and standard deviation
of Spearman’s $\rho$ across all Java methods where correlation is
statistically significant ($p\leq 0.05$), and (2) the proportion of Java
methods demonstrating a statistically significant correlation ($p\leq 0.05$).

In addition, we group all AST tokens into 18 semantic categories (e.g., method
call, operator, etc.) and investigate how much humans333We use fixation
durations to represent human focus. We empirically verify that using fixation
count yields similar results. and language models focus on each semantic
category. The focus score assigned to each semantic category is the sum of the
focus scores assigned to each AST token belonging to that semantic category.
To counter biases where certain semantic categories contain more AST tokens or
appear more frequently, we report the relative difference between machine and
human foci for each semantic category. That is, we average the six language
models’ focus scores per category and report
$|\frac{foc_{machine}-foc_{human}}{foc_{human}}|$.

####  3.3.2. RQ2

Here, a human expert provides quality ratings for summaries generated by each
language model for every Java method using four criteria: accuracy,
completeness, conciseness, and readability. Next, we calculate the Spearman’s
$\rho$ between each language model’s focus vector and the human fixation
duration vector across all Java methods where correlation is significant
($p\leq 0.05$). We then append all such statistically significant $\rho$’s to
form a vector, denoted $v_{cor}$, to represent the degrees of alignment
between machine and human foci across the Java methods investigated444Note
that $v_{cor}$ contains Spearman’s $\rho$’s obtained from all six language
models. We empirically verify that conducting the analogous analysis for each
language model separately yields a similar result..

Subsequently, we determine whether this alignment is correlated with the rated
quality of summaries. Specifically, we construct four other vectors,
$\\{v_{acc},v_{com},v_{con},v_{rea}\\}$, containing the accuracy,
completeness, conciseness, and readability scores respectively. At each index
$i$, $\\{v_{cor}[i],v_{acc}[i],v_{com}[i],v_{con}[i],v_{rea}[i]\\}$ are
respectively the Spearman’s $\rho$, summary accuracy, completeness,
conciseness, and readability of the same language model applied on the same
Java method. We then measure and report the Spearman’s rank correlation
between $v_{cor}$ and $v_{i}$, where
$v_{i}\in\\{v_{acc},v_{com},v_{con},v_{rea}\\}$.

##  4\. Results

###  4.1. RQ1: General Correlation

As shown in Table 1, there is a general lack of correlation between human and
machine foci. We highlight that the means and standard deviations in Table 1
are only calculated from Java methods where the Spearman’s $\rho$ is
statistically significant (with $p\leq 0.05$). In practice, between any pair
of human-LLM focus sources, at most $22\%$ of the 68 Java methods yield a
Spearman’s $\rho$ with $p\leq 0.05$. As a baseline, the Spearman’s $\rho$ has
$p\leq 0.05$ for all Java methods between human duration and fixation focus
vectors, and for $85\%$ of Java methods between any two language model’s focus
vectors. This implies that any existing correlation between human and machine
foci is not widespread across the Java methods studied.

Furthermore, among those Java methods where the correlation is statistically
significant, the mean Spearman’s $\rho$ is small with a large standard
deviation. In fact, for most such methods where a model and human show
significant correlation in focus, the Spearman’s $\rho$ is often either around
$0.5$ or $-0.5$, but rarely in between. This suggests the relationship between
human and machine foci varies significantly depending on the specific Java
method.

Interestingly, although few-shot-alignment in GPT-few-shot renders the model’s
generated summaries more similar to those of humans, this does not lead to
higher correlations between model and human foci. In addition, feature
attribution in all language models is moderately or strongly positively
correlated with each other on a majority of Java methods, which intuitively
makes sense since all six models are based on the Transformer architecture.

We also investigate how language models’ focus on each semantic category
differs from that of humans. As shown in Figure 1, language models’ generation
of code summaries seems to be more reliant on comments, return values, and
specific statements such as literals and assignments, and less reliant on
method calls and variables/methods not defined explicitly within the Java
method.

Discussion Point 1: We find no evidence that feature attribution in language
models is correlated with programmers’ visual focus. Several possible
interpretations can be inferred: (1) Alternative methods may be needed to
assess feature influence in black-box language models for code summarization,
aiming for better alignment with human attention. (2) Access to the internal
workings of proprietary models might become critical if white-box models offer
more human-aligned insights into explainable language models for code
(paltenghi2022extracting, ). (3) It is possible that language models and
humans reason about code differently when summarizing source code.

Table 1. Pair-wise correlation among focus sources; “Duration” and “Count” are human visual focus. Each cell shows the means and standard deviations of Spearman’s $\rho$ for all Java methods showing significant correlation ($p\leq 0.05$). | Duration | Count | GPT4 | GPT-few | GPT3.5 | StarCoder | CodeLlama | NCS
---|---|---|---|---|---|---|---|---
Duration | 1.00$\pm$0.00 | 0.88$\pm$0.06 | -0.11$\pm$0.41 | -0.13$\pm$0.42 | -0.09$\pm$0.52 | -0.18$\pm$0.48 | -0.18$\pm$0.42 | -0.24$\pm$0.40
Count | — | 1.00$\pm$0.00 | 0.01$\pm$0.45 | -0.24$\pm$0.33 | -0.10$\pm$0.48 | -0.31$\pm$0.29 | -0.13$\pm$0.43 | -0.33$\pm$0.33
GPT4 | — | — | 1.00$\pm$0.00 | 0.68$\pm$0.12 | 0.76$\pm$0.12 | 0.67$\pm$0.14 | 0.67$\pm$0.14 | 0.55$\pm$0.13
GPT-few | — | — | — | 1.00$\pm$0.00 | 0.72$\pm$0.12 | 0.62$\pm$0.15 | 0.64$\pm$0.15 | 0.55$\pm$0.13
GPT3.5 | — | — | — | — | 1.00$\pm$0.00 | 0.65$\pm$0.16 | 0.67$\pm$0.15 | 0.58$\pm$0.13
StarCoder | — | — | — | — | — | 1.00$\pm$0.00 | 0.66$\pm$0.15 | 0.59$\pm$0.11
Code Llama | — | — | — | — | — | — | 1.00$\pm$0.00 | 0.56$\pm$0.14
NCS | — | — | — | — | — | — | — | 1.00$\pm$0.00
Figure 1. How much more/less do language models focus on each semantic
category compared to humans?  Figure 2. Average ratings of model-generated
summaries. Table 2. Correlation between human-machine focus alignment and
summary quality (assessed by four metrics).

| Accuracy | Completeness | Conciseness | Readability
---|---|---|---|---
Spearman’s $\rho$ | -0.1279 | 0.1309 | 0.0194 | -0.0717
$p$-value | 0.3862 | 0.3753 | 0.8960 | 0.6280

###  4.2. RQ2: Summary Qualities

There is also a lack of correlation between the quality of summaries generated
by language models and how well their focus on code aligns with humans’. The
large p-values in Table 1 suggest that, regardless of which metric is used to
assess summary quality, there is a lack of statistically significant
correlation between the quality of a model-generated summary on a Java method
and how well the model’s focus aligns with that of humans on that Java method.
Furthermore, Figure 2 shows that NCS produces worse summaries than the other
five models. Although Table 1 seems to suggest that NCS’s focus is more
negatively aligned with human attention, we find no statistically significant
metrics supporting such a claim, partially due to the small sample size of
Java methods yielding statistically significant Spearman’s $\rho$.

In general, Table 1 suggests that feature attribution in NCS is still
moderately positively aligned with that in other language models on a majority
of Java methods. This indicates the likelihood that aspects other than feature
attribution are more indicative of and critical to a language model’s
performance in code summarization.

Discussion Point 2: With a substantial body of work in NLP showing that
aligning neural models with human visual patterns can lead to performance
improvement (sood2020improving, ; zhang2019using, ; klerke2016improving, ;
barrett2018sequence, ), we contain our conclusion to the SHAP measure of
feature attribution and the human attention as measured in an eye-tracking
experiment. The link between human attention and feature attribution to
machine models is a subject of intense scientific investigation. We contribute
to the debate with this finding that SHAP did not correlate with human eye
attention in the measures or models we studied.

##  5\. Conclusion

In this paper, we use a state-of-the-art, black-box, perturbation-based method
to assess feature attribution in language models on code summarization tasks.
We then compare the model-determined important AST tokens with those
identified by human visual focus, as measured through eye-tracking. The
results suggest that using SHAP to measure feature attribution does not
provide explainability of language models through establishing correlations
between machine and human foci. Generally, our work can be interpreted in two
ways. First, feature attribution measured by SHAP may not be the best way to
interpret a language model’s focus during code summarization as it fails to
establish similarities with human focus. Alternatively, it may be the case
that machines reason about code differently from humans when tasked to
summarize source code.

##  6\. Acknowledgement

This research was supported by NSF CCF-2211429, NSF CCF-2211428, NSF
CCF-2100035 and NSF SaTC-2312057.

## References

  * (1) Tobii pro fusion user manual, Jun 2023.
  * (2) Abid, N. J., Maletic, J. I., and Sharif, B. Using developer eye movements to externalize the mental model used in code summarization tasks.  In Proceedings of the 11th ACM Symposium on Eye Tracking Research & Applications (2019), pp. 1–9.
  * (3) Ahmad, W. U., Chakraborty, S., Ray, B., and Chang, K.-W. A transformer-based approach for source code summarization.  arXiv preprint arXiv:2005.00653 (2020).
  * (4) Bansal, A., Su, C.-Y., Karas, Z., Zhang, Y., Huang, Y., Li, T. J.-J., and McMillan, C. Modeling programmer attention as scanpath prediction.  arXiv preprint arXiv:2308.13920 (2023).
  * (5) Barrett, M., Bingel, J., Hollenstein, N., Rei, M., and Søgaard, A. Sequence classification with human attention.  In Proceedings of the 22nd conference on computational natural language learning (2018), pp. 302–312.
  * (6) Collard, M. L., Decker, M. J., and Maletic, J. I. srcml: An infrastructure for the exploration, analysis, and manipulation of source code: A tool demonstration.  In 2013 IEEE International conference on software maintenance (2013), IEEE, pp. 516–519.
  * (7) Galassi, A., Lippi, M., and Torroni, P. Attention in natural language processing.  IEEE transactions on neural networks and learning systems 32, 10 (2020), 4291–4308.
  * (8) Hooker, S., Erhan, D., Kindermans, P.-J., and Kim, B. Evaluating feature importance estimates.
  * (9) Huber, D., Paltenghi, M., and Pradel, M. Where to look when repairing code? comparing the attention of neural models and developers.  arXiv preprint arXiv:2305.07287 (2023).
  * (10) Klerke, S., Goldberg, Y., and Søgaard, A. Improving sentence compression by learning to predict gaze.  arXiv preprint arXiv:1604.03357 (2016).
  * (11) Kou, B., Chen, S., Wang, Z., Ma, L., and Zhang, T. Is model attention aligned with human attention? an empirical study on large language models for code generation.  arXiv preprint arXiv:2306.01220 (2023).
  * (12) LeClair, A., and McMillan, C. Recommendations for datasets for source code summarization.  arXiv preprint arXiv:1904.02660 (2019).
  * (13) Li, R., Allal, L. B., Zi, Y., Muennighoff, N., Kocetkov, D., Mou, C., Marone, M., Akiki, C., Li, J., Chim, J., et al. Starcoder: may the source be with you!  arXiv preprint arXiv:2305.06161 (2023).
  * (14) Liu, Y., Tantithamthavorn, C., Liu, Y., and Li, L. On the reliability and explainability of automated code generation approaches.  arXiv preprint arXiv:2302.09587 (2023).
  * (15) Lundberg, S. M., and Lee, S.-I. A unified approach to interpreting model predictions.  Advances in neural information processing systems 30 (2017).
  * (16) Molnar, C. Interpretable machine learning.  Lulu. com, 2020.
  * (17) Nijkamp, E., Pang, B., Hayashi, H., Tu, L., Wang, H., Zhou, Y., Savarese, S., and Xiong, C. Codegen: An open large language model for code with multi-turn program synthesis.  arXiv preprint arXiv:2203.13474 (2022).
  * (18) OpenAI.  Introducing chatgpt.  https://openai.com/blog/chatgpt/, 2023.  Accessed: 11/20/2023.
  * (19) OpenAI, R. Gpt-4 technical report. arxiv 2303.08774.  View in Article 2 (2023).
  * (20) Paltenghi, M., Pandita, R., Henley, A. Z., and Ziegler, A. Extracting meaningful attention on source code: An empirical study of developer and neural model code exploration.  arXiv preprint arXiv:2210.05506 (2022).
  * (21) Paltenghi, M., and Pradel, M. Thinking like a developer? comparing the attention of humans with neural models of code.  In 2021 36th IEEE/ACM International Conference on Automated Software Engineering (ASE) (2021), IEEE, pp. 867–879.
  * (22) Rodeghero, P., and McMillan, C. An empirical study on the patterns of eye movement during summarization tasks.  In 2015 ACM/IEEE International Symposium on Empirical Software Engineering and Measurement (ESEM) (2015), IEEE, pp. 1–10.
  * (23) Roziere, B., Gehring, J., Gloeckle, F., Sootla, S., Gat, I., Tan, X. E., Adi, Y., Liu, J., Remez, T., Rapin, J., et al. Code llama: Open foundation models for code.  arXiv preprint arXiv:2308.12950 (2023).
  * (24) Sharafi, Z., Shaffer, T., Sharif, B., and Guéhéneuc, Y.-G. Eye-tracking metrics in software engineering.  In 2015 Asia-Pacific Software Engineering Conference (APSEC) (2015), IEEE, pp. 96–103.
  * (25) Shrikumar, A., Greenside, P., and Kundaje, A. Learning important features through propagating activation differences.  In International conference on machine learning (2017), PMLR, pp. 3145–3153.
  * (26) Sood, E., Tannert, S., Müller, P., and Bulling, A. Improving natural language processing tasks with human gaze-guided neural attention.  Advances in Neural Information Processing Systems 33 (2020), 6327–6341.
  * (27) Spearman, C. The proof and measurement of association between two things.
  * (28) Sundararajan, M., Taly, A., and Yan, Q. Axiomatic attribution for deep networks.  In International conference on machine learning (2017), PMLR, pp. 3319–3328.
  * (29) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I. Attention is all you need.  Advances in neural information processing systems 30 (2017).
  * (30) Wang, Y., Wang, K., and Wang, L. Wheacha: A method for explaining the predictions of code summarization models.  arXiv preprint arXiv:2102.04625 (2021).
  * (31) Xu, F. F., Alon, U., Neubig, G., and Hellendoorn, V. J. A systematic evaluation of large language models of code.  In Proceedings of the 6th ACM SIGPLAN International Symposium on Machine Programming (2022), pp. 1–10.
  * (32) Zeng, Z., Tan, H., Zhang, H., Li, J., Zhang, Y., and Zhang, L. An extensive study on pre-trained models for program understanding and generation.  In Proceedings of the 31st ACM SIGSOFT international symposium on software testing and analysis (2022), pp. 39–51.
  * (33) Zhang, K., Li, G., and Jin, Z. What does transformer learn about source code?  arXiv preprint arXiv:2207.08466 (2022).
  * (34) Zhang, Y., and Zhang, C. Using human attention to extract keyphrase from microblog post.  In Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics (2019), pp. 5867–5872.
  * (35) Zhang, Z., Zhang, H., Shen, B., and Gu, X. Diet code is healthy: Simplifying programs for pre-trained models of code.  In Proceedings of the 30th ACM Joint European Software Engineering Conference and Symposium on the Foundations of Software Engineering (2022), pp. 1073–1084.