# Bispectrum-based Cross-frequency Functional Connectivity: Classification of
Alzheimer’s disease

Dominik Klepl  Centre for Computational Science and Mathematical Modelling,
Coventry University, Coventry CV1 2JH, UK  Fei He<EMAIL_ADDRESS><EMAIL_ADDRESS>Centre for Computational Science and Mathematical
Modelling, Coventry University, Coventry CV1 2JH, UK  Min Wu  Institute for
Infocomm Research, Agency for Science, Technology and Research (A*STAR),
138632, Singapore  Daniel J. Blackburn  Department of Neuroscience, University
of Sheffield, Sheffield, S10 2HQ, UK  Ptolemaios G. Sarrigiannis  Department
of Neurophysiology, Royal Devon and Exeter NHS Foundation Trust, Exeter, EX2
5DW, UK

###### Abstract

Alzheimer’s disease (AD) is a neurodegenerative disease known to affect brain
functional connectivity (FC). Linear FC measures have been applied to study
the differences in AD by splitting neurophysiological signals such as
electroencephalography (EEG) recordings into discrete frequency bands and
analysing them in isolation. We address this limitation by quantifying cross-
frequency FC in addition to the traditional within-band approach. Cross-
bispectrum, a higher-order spectral analysis, is used to measure the nonlinear
FC and is compared with the cross-spectrum, which only measures the linear FC
within bands. Each frequency coupling is then used to construct an FC network,
which is in turn vectorised and used to train a classifier. We show that
fusing features from networks improves classification accuracy. Although both
within-frequency and cross-frequency networks can be used to predict AD with
high accuracy, our results show that bispectrum-based FC outperforms cross-
spectrum suggesting an important role of cross-frequency FC.

Clinical relevance— This establishes diagnostic relevance of cross-frequency
coupling in Alzheimer’s disease.

##  I Introduction

Alzheimer’s disease (AD) causes early neural degradation leading to cell death
and synaptic loss, and is the most common form of dementia. Studies have shown
that AD alters functional connectivity (FC) over multiple scales
pievani2011functional ; konig2005decreased ; jeong2004eeg .
Electroencephalography (EEG) is a commonly used method to study and classify
AD. The main EEG characteristics associated with AD are slowing of signals and
decrease in synchronisation jeong2004eeg ; ghorbanian2015slowing ;
konig2005decreased ; babiloni2016synchronisationreview . Slowing of EEG in AD
was observed as an increased activity in $\delta$ and $\theta$ frequency bands
and decreased activity in $\alpha$, $\beta$ and $\gamma$ frequency bands
jeong2004eeg ; ghorbanian2015slowing . Similarly, AD shows changes in
synchronisation within low-frequency bands ($<$12 Hz) and is associated with
altered FC babiloni2016synchronisationreview . However, these characteristics
are typically measured on a single EEG or pairwise channel level. In contrast,
network based FC analysis can reveal additional characteristics of AD across
multiple channels and scales kabbara2018integrationsegregation ;
dai2019disrupted . However, these studies focus only on specific frequency
bands, i.e. the within-frequency coupling (WFC). WFC networks of AD have been
analysed using various FC methods such as coherence
blinowska2017classification_coherence , wavelet coherence
jeong2016waveletcoherence  and phase lag index (PLI)
nobukawa2020classification_PLI . WFC networks were also used for
classification of AD using a multivariate autoregressive model
blinowska2017classification_coherence  and support vector machine (SVM)
nobukawa2020classification_PLI . ††©2022 IEEE. Personal use of this material
is permitted. Permission from IEEE must be obtained for all other uses, in any
current or future media, including reprinting/republishing this material for
advertising or promotional purposes, creating new collective works, for resale
or redistribution to servers or lists, or reuse of any copyrighted component
of this work in other works.

Figure 1: A conceptual schematic of constructing the proposed cross-bispectrum
(CBS) networks. First, each EEG signal is cleaned and scaled (left). For each
pair of EEG electrodes, CBS is estimated (middle). The frequency bands
coupling edge weights are given by the maximum value within the respective CBS
window. From each CBS, 25 edges are inferred. This way, 5 within-frequency
(WFC) and 20 cross-frequency (CFC) networks are constructed (right).

This study aims to extend the WFC FC analysis by taking the cross-frequency
coupling (CFC) jirsa2013crossbispectrumcoupling  into account and investigate
the predictive power of CFC for AD diagnosis. Only one CFC measure, i.e. phase
synchronisation index (PSI), had been used to reconstruct CFC networks of AD
cai2018ADcrossphase . However, this study did not use any classification
methods. Moreover, PSI measures only phase-phase CFC.

(Cross-) Bispectrum (CBS), a higher-order spectral analysis, quantifies
quadratic coupling between two frequency components and their algebraic sum
HE2021 . It has been shown to detect phase-amplitude and phase-phase coupling
jirsa2013crossbispectrumcoupling . CBS-based EEG features have been proposed
as biomarkers of epilepsy mahmoodian2019epilepsybispectrum , Parkinson’s
disease yuvaraj2018PDindexbispectra  and AD wang2015multiple . However, these
applications study only a few EEG channels. In contrast, our study computes
CBS between all pairs of EEG channels to construct both WFC and CFC FC
networks (Fig. 1).

The CBS estimates of FC are computed to investigate the differences between AD
and healthy controls (HC) in nonlinear WFC and CFC, compared to linear WFC
measured with cross-spectrum (CS). To the best of our knowledge, this is the
first application of CBS in reconstructing FC networks of AD. Finally, we use
the estimated networks to train classifiers, in order to examine the
predictive power of WFC and CFC.

##  II Data

EEG recordings were collected from 20 AD patients and 20 HC under 70. A
detailed description of the experimental design and confirmation of the
diagnosis is provided in blackburn2018synchronisation . AD patients that were
in the mild to moderate stage with the average Mini Mental State Examination
score of 20.1 $(sd=4)$ were recruited in the Sheffield Teaching Hospital. Age
and gender-matched HC participants were recruited. A consent was obtained from
all participants.

EEG was acquired using an XLTEK 128-channel headbox, Ag/AgCL electrodes with a
sampling frequency of 2 kHz using a modified 10-10 overlapping and 10-20
international electrode placement system with a referential montage with a
linked earlobe reference blackburn2018synchronisation . 30 minutes long
sessions at rest were recorded, consisting of alternating two-minute-long
epochs during which the participants had their eyes closed (EC) or open. Only
the EC epochs are used in this study.

All the recordings were reviewed by a neurophysiologist on the XLTEK review
station with time-locked video recordings. Three 12-second-long artefact-free
epochs per subject were isolated. Finally, 23 bipolar channels were created:
F8–F4, F7–F3, F4–C4, F3–C3, F4–FZ, FZ–CZ, F3–FZ, T4–C4, T3–C3, C4–CZ, C3–CZ,
CZ–PZ, C4–P4, C3–P3, T4–T6, T3–T5, P4–PZ, P3–PZ, T6–O2, T5–O1, P4–O2, P3–O1
and O1–O2 blackburn2018synchronisation .

###  II.1 EEG pre-processing

EEG signals were band-pass filtered to 0.1-100 Hz range using a zero-phase 5th
order Butterworth filter. 50 Hz power line noise was removed using a zero-
phase 4th order Butterworth stop-band filter; and the data were down-sampled
to 250 Hz using an 8th order Chebyshev type I filter. The signals were
normalised to zero mean and unit standard deviation.

Figure 2: Average connectivity matrices measured with cross-spectrum of (A) AD
and (B) HC. For visualisation, the values were min-max normalised.  Figure 3:
Average connectivity matrices measured with cross-bispectrum of (A) AD and (B)
HC with input frequency on the vertical facets and output frequency on
horizontal. For visualisation purposes, the values were min-max normalised.

##  III Methods

###  III.1 Cross-spectrum and cross-bispectrum

The CS of signal $S_{X}$ is calculated via smoothed periodogram. Fast Fourier
Transform (FFT) is used to estimate the periodogram with Daniell smoothers.
The periodogram is computed over 256 frequency bins (0.98 Hz bandwidth). CS at
frequency $f$ is then computed as:

$CS_{XY}(f)=S_{X}(f)\cdot S_{Y}(f).$ (1)

A direct FFT-based method is used to estimate CBS:

$CBS_{XY}(f_{1},f_{2})=\langle X_{T}(f_{1})\cdot Y_{T}(f_{2})\cdot
Y_{T}^{\ast}(f_{1}+f_{2})\rangle,$ (2)

where $\langle\cdot\rangle$ denotes averaging, $X_{T}(f)$ is a Fourier
Transform of signal $X$ over an interval $T$ and $Y^{\ast}$ is the complex
conjugate. 256-point FFT is used. CBS is computed over 50 segments with 50%
overlap. The estimated CBS is then smoothed in frequency using a Rao-Gabr
window (size 5).

CS and CBS were computed for all pairs of EEG channels. Five frequency bands
$b$ are considered: $\delta$ $(0.5-5Hz)$, $\theta$ $(5-8Hz)$, $\alpha$
$(8-16Hz)$, $\beta$ $(16-32Hz)$ and $\gamma$ $(32-100Hz)$.

The connectivity $CN$ between channels $X$ and $Y$ and frequency bands $b_{X}$
and $b_{Y}$ is computed as:

$CN_{XY}^{CS}(b_{X})=\max(|CS_{XY}(f\in b_{X})|),b_{X}=b_{Y},$ (3)
$CN_{XY}^{CBS}(b_{X},b_{Y})=\max(|CBS_{XY}(f_{1}\in b_{X},f_{2}\in b_{Y})|),$
(4)

for CS and CBS respectively. This results in 5 WFC and 20 CFC measures for
each pair of channels. Note that CBS edges are directed, thus we collapse them
to undirected by taking the average of edge weights.

To reduce the risk of false positives induced by spurious random coupling,
surrogate thresholding was used theiler1992surrogate . For each pair of
channels, 200 surrogate signals were generated using the Fourier transform
surrogate method, which scrambles the phase, and their CS and CBS are
computed. The 95% confidence interval of surrogate values is used as a
threshold. Coupling values below the threshold are set to zero. We obtain a
set of brain networks for each EEG epoch, i.e. $N\times N$ matrices ($N=23$).
For CS and CBS, there are 5 and 25 connectivity matrices, respectively.

To further reduce the risk of false-negative edges, the weighted networks are
filtered using relative quantile-based thresholding which preserves only the
top-k % of the strongest edges ($k\in\\{10,20,30,40,100\\}$.

###  III.2 Network classification

Finally, we use the networks to elucidate their discriminatory power for
classification of AD. The upper triangular of adjacency matrix of each network
is vectorised and forms a feature vector used to train a classifier. SVM with
radial basis kernel is used and principal component analysis is employed, such
that the components capture 90% variance of the data. Moreover, the features
are scaled to zero mean and unit standard deviation. A leave-one-subject-out
cross-validation is used in this work.

Two categories of classifiers are trained: network level and concatenated-
networks level. Network level classifiers are fitted for each constructed
network and each network filter to test which frequency coupling types can be
used for classification, i.e. show the largest differences between AD and HC.
Thus, we fit 125 and 25 (frequency band combinations $\times$ network filters)
classifiers using CBS and CS respectively.

Concatenated-networks level classifiers are fitted by concatenating the
vectorised adjacency matrices for each network filter separately. Thus, the
feature vectors consists of concatenated adjacency matrices of 25 and 5
networks for CS and CBS, respectively. We hypothesise that the information
captured by each network is at least partially unique, thus a classifier
trained on the concatenated networks should outperform the classifiers trained
on individual networks, as it can leverage the information from all networks.

##  IV Results and Discussion

Figure 4: Accuracy of classifiers trained across frequency couplings and
filtered networks using cross-spectrum (orange) and cross-bispectrum (blue)
networks.

###  IV.1 Connectivity matrices

Averaged connectivity matrices (Fig. 2 and 3) indicate that both CS and CBS
seem to detect the differences in the networks. In CS, these differences are
most pronounced in $\delta$ and $\theta$ bands where AD have higher FC, and in
$\beta$ where AD have lower FC. In CBS, differences can be observed in almost
all frequency bands and their couplings. The differences are rather obvious in
$\delta$, $\theta$, $\alpha$ and their CFC. Differences in these bands (both
WFC and CFC) were reported by Cai et al. cai2018ADcrossphase  using PSI, but
in the opposite direction to our findings. This might be explained by the
differences between CBS and PSI, as PSI is a purely phase measure unlike the
CBS jirsa2013crossbispectrumcoupling . Our findings are consistent with the
literature reporting the increased activity in $\delta$ and $\theta$ in AD
jeong2004eeg ; konig2005decreased . Also, the distortion of structure observed
in multiple CS and CBS networks is consistent with the disrupted information
processing in AD dai2019disrupted .

###  IV.2 Classification results

First, SVM classifiers were trained for each network separately to compare the
discriminatory power of CS and CBS based networks (Fig. 4). However, a direct
comparison of CS and CBS is only possible for the WFC classifiers (along the
diagonal of Fig. 4). CS-based models perform better in $\delta$ and $\theta$.
On the other hand, CBS-based models perform better in $\alpha$, $\beta$ and in
$\gamma$ (100%). Thus, in the WFC comparison, neither of the methods performs
consistently better, suggesting that CBS-based networks are equally capable of
capturing the differences in WFC as the CS-based networks.

Additionally, the results show that multiple types of CBS-based CFC networks
can deliver comparable or even slightly better accuracy (off diagonal of Fig.
4). This is most clear in the $\delta-\theta$, $\delta-\alpha$,
$\theta-\delta$, $\theta-\alpha$, $\alpha-\delta$ and $\alpha-\theta$ models.
This highlights the importance of studying CFC networks in addition to the WFC
networks. Next, we trained concatenated-networks models to examine whether
leveraging information from multiple networks improves the classification
performance. This is confirmed in Table 1. All CBS-based models outperform the
CS-based models suggesting that including information about nonlinear and CFC
coupling can be crucial for modelling and classification of AD.

Table 1: Performance of concatenated-networks models Model | Accuracy | Sensitivity | Specificity
---|---|---|---
Spectrum - Top 100% | 73% | 71.8% | 74.2%
Spectrum - Top 40% | 72.4% | 73.7% | 71.4%
Spectrum - Top 30% | 72.4% | 74.3% | 71%
Spectrum - Top 20% | 71.3% | 72.4% | 70.4%
Spectrum - Top 10% | 70.7% | 80% | 66.4%
Bispectrum - Top 100% | 74.1% | 72.9% | 75.3%
Bispectrum - Top 40% | 73.6% | 72.6% | 74.4%
Bispectrum - Top 30% | 75.9% | 72.3% | 80%
Bispectrum - Top 20% | 77.6% | 74.2% | 81.5%
Bispectrum - Top 10% | 77% | 81.4% | 74%

###  IV.3 Conclusions and Future work

We have demonstrated that CBS and CS detect similar differences between AD and
HC networks, but CBS has an advantage of capturing cross-frequency and
nonlinear interactions. Although CBS was shown to be a powerful tool to detect
various types of WFC and CFC, such as phase-phase and phase-amplitude, it is
not possible to distinguish between them. Therefore, a combination of CBS with
other types of CFC methods could be a plausible direction for future research.

Finally, we use relatively simple features for the classification purpose. The
use of feature selection methods might further improve the classification
performance by removing uninformative and correlated features. In future
research, it might also be important and interesting to explore graph-based
features that would capture the differences between AD and HC in a lower-
dimensional space more efficiently.

## Acknowledgement

The original data collection was funded by a grant from the Alzheimer’s
Research UK (ARUK-PPG20114B-25). The views expressed are those of the
author(s) and not necessarily those of the NHS, the NIHR or the Department of
Health.

## References

  * [1] Michela Pievani, Willem de Haan, Tao Wu, William W Seeley, and Giovanni B Frisoni.  Functional network disruption in the degenerative dementias.  The Lancet Neurol., 10(9):829–843, 2011.
  * [2] Thomas König, L Prichep, Thomas Dierks, Daniela Hubl, LO Wahlund, ER John, and V Jelic.  Decreased EEG synchronization in Alzheimer’s disease and mild cognitive impairment.  Neurobiol. Aging, 26(2):165–171, 2005.
  * [3] Jaeseung Jeong.  EEG dynamics in patients with Alzheimer’s disease.  Clin. Neurophysiol., 115(7):1490–1505, 2004.
  * [4] Parham Ghorbanian, David M Devilbiss, Terry Hess, Allan Bernstein, Adam J Simon, and Hashem Ashrafiuon.  Exploration of EEG features of Alzheimer’s disease using continuous wavelet transform.  Med. Biol. Eng., 53(9):843–855, 2015.
  * [5] Claudio Babiloni, Roberta Lizio, Nicola Marzano, Paolo Capotosto, Andrea Soricelli, Antonio Ivano Triggiani, Susanna Cordone, Loreto Gesualdo, and Claudio Del Percio.  Brain neural synchronization and functional coupling in Alzheimer’s disease as revealed by resting state EEG rhythms.  Int. J. Psychophysiol., 103:88–102, 2016.
  * [6] A Kabbara, H Eid, W El Falou, M Khalil, F Wendling, and M Hassan.  Reduced integration and improved segregation of functional brain networks in Alzheimer’s disease.  J. Neural. Eng., 15(2):026023, 2018.
  * [7] Zhengjia Dai, Qixiang Lin, Tao Li, Xiao Wang, Huishu Yuan, Xin Yu, Yong He, and Huali Wang.  Disrupted structural and functional brain networks in Alzheimer’s disease.  Neurobiol. Aging, 75:71–82, 2019.
  * [8] Katarzyna J Blinowska, Franciszek Rakowski, Maciej Kaminski, Fabrizio De Vico Fallani, Claudio Del Percio, Roberta Lizio, and Claudio Babiloni.  Functional and effective brain connectivity for discrimination between Alzheimer’s patients and healthy individuals: A study on resting state EEG rhythms.  Clin. Neurophysiol., 128(4):667–680, 2017.
  * [9] Dong-Hwa Jeong, Young-Do Kim, In-Uk Song, Yong-An Chung, and Jaeseung Jeong.  Wavelet energy and wavelet coherence as EEG biomarkers for the diagnosis of Parkinson’s disease-related dementia and Alzheimer’s disease.  Entropy, 18(1):8, 2016.
  * [10] Sou Nobukawa, Teruya Yamanishi, Shinya Kasakawa, Haruhiko Nishimura, Mitsuru Kikuchi, and Tetsuya Takahashi.  Classification methods based on complexity and synchronization of electroencephalography signals in Alzheimer’s disease.  Front. Psychiatry, 11:255, 2020.
  * [11] Viktor Jirsa and Viktor Müller.  Cross-frequency coupling in real and virtual brain networks.  Front. Computat. Neurosci., 7:78, 2013.
  * [12] Lihui Cai, Xile Wei, Jiang Wang, Haitao Yu, Bin Deng, and Ruofan Wang.  Reconstruction of functional brain network in Alzheimer’s disease via cross-frequency phase synchronization.  Neurocomputing, 314:490–500, 2018.
  * [13] Fei He and Yuan Yang.  Nonlinear system identification of neural systems from neurophysiological signals.  Neuroscience, 458:213–228, 2021.
  * [14] Naghmeh Mahmoodian, Axel Boese, Michael Friebe, and Javad Haddadnia.  Epileptic seizure detection using cross-bispectrum of electroencephalogram signal.  Seizure, 66:4–11, 2019.
  * [15] Rajamanickam Yuvaraj, U Rajendra Acharya, and Yuki Hagiwara.  A novel Parkinson’s disease diagnosis index using higher-order spectra features in EEG signals.  Neural Comput. Appl., 30(4):1225–1235, 2018.
  * [16] Ruofan Wang, Jiang Wang, Shunan Li, Haitao Yu, Bin Deng, and Xile Wei.  Multiple feature extraction and classification of electroencephalograph signal for Alzheimers’ with spectrum and bispectrum.  Chaos: An Interdisciplinary J. of Nonlinear Science, 25(1):013110, 2015.
  * [17] Daniel J. Blackburn, Yifan Zhao, Matteo De Marco, Simon M. Bell, Fei He, Hua-Liang Wei, Sarah Lawrence, Zoe C. Unwin, Michelle Blyth, and Jenna Angel.  A pilot study investigating a novel non-linear measure of eyes open versus eyes closed EEG synchronization in people with Alzheimer’s disease and healthy controls.  Brain Sci., 8(7):134, 2018.
  * [18] James Theiler, Stephen Eubank, André Longtin, Bryan Galdrikian, and J Doyne Farmer.  Testing for nonlinearity in time series: the method of surrogate data.  Physica D: Nonlinear Phenomena, 58(1-4):77–94, 1992.


# Swimming, Feeding and Inversion of Multicellular Choanoflagellate Sheets

Lloyd Fung<EMAIL_ADDRESS>Department of Applied Mathematics and Theoretical
Physics, Centre for Mathematical Sciences,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom
Adam Konkol<EMAIL_ADDRESS>Department of Applied Mathematics and
Theoretical Physics, Centre for Mathematical Sciences,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom
Takuji Ishikawa<EMAIL_ADDRESS>Department of Biomedical Engineering,
Tohoku University, 6-6-01 Aoba, Aramaki, Aoba-ku, Sendai 980-8579, Japan  Ben
T. Larson<EMAIL_ADDRESS>Department of Biochemistry & Biophysics,
University of California, San Francisco, 600 16th St., San Francisco, CA
94143-2200, USA
Thibaut Brunet<EMAIL_ADDRESS>Department of Cell Biology and
Infection, Institut Pasteur, 25-28 rue du Dr. Roux, 75724 Paris Cedex 15,
France  Raymond E. Goldstein<EMAIL_ADDRESS>Department of
Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom

###### Abstract

The recent discovery of the striking sheet-like multicellular choanoflagellate
species Choanoeca flexa that dynamically interconverts between two
hemispherical forms of opposite orientation raises fundamental questions in
cell and evolutionary biology, as choanoflagellates are the closest living
relatives of animals. It similarly motivates questions in fluid and solid
mechanics concerning the differential swimming speeds in the two states and
the mechanism of curvature inversion triggered by changes in the geometry of
microvilli emanating from each cell. Here we develop fluid dynamical and
mechanical models to address these observations and show that they capture the
main features of the swimming, feeding, and inversion of C. flexa colonies.

Some of the most fascinating processes in the developmental biology of complex
multicellular organisms involve radical changes in geometry or topology. From
the folding of tissues during gastrulation gastrulation  to the formation of
hollow spaces in plants Goriely , these processes generally involve
coordinated cell shape changes, cellular division, migration and apoptosis,
and formation of an extracellular matrix (ECM). It has become clear through
multiple strands of research that evolutionary precedents for these processes
exist in some of the simplest multicellular organisms such as green algae
hohn2016distinct ; Volvox_prl  and choanoflagellates BrunetKing , the latter
being the closest living relatives of animals. Named for their funnel-shaped
collar of microvilli that facilitates filter feeding from the flows driven by
their beating flagellum, choanoflagellates serve as model organisms for the
study of the evolution of multicellularity.

While well-known multicellular choanoflagellates exist as linear chains or
“rosettes” fairclough2010multicellular  held together by an ECM
larson2020biophysical , a new species named Choanoeca flexa was recently
discovered brunet2019light  with an unusual sheet-like geometry (Fig. 1) in
which hundreds of cells adhere to each other by the tips of their microvilli,
without an ECM leadbeater1983life . The sheets can exist in two forms with
opposite curvature, one with flagella pointing towards the center of curvature
[“flag-in”] with a relatively large spacing between cells, and another with
the opposite arrangement [“flag-out”] with more tightly-packed cells.
Transformations from flag-in to flag-out can be triggered by darkness, and
occur in $\sim\\!10\,$s. Compared to the flag-out form, the flag-in state has
limited motility and is better suited to filter-feeding. It was conjectured
brunet2019light  that the darkness-induced transition to the more motile form
is a type of photokinesis.

Figure 1: The multicellular choanoflagellate Choanoeca flexa. Top views of (a)
flag-in and (b) flage-out states at times relative to removal of light. (c)
Close-up of the collar connections in the two states. (d) Electron micrograph
showing round white cell bodies connected by microvilli. Adapted from
brunet2019light . Figure 2: Models for _C. flexa_. (a) Cell body and flagellar
force in the flag-in state. (b,c) Mechanical model of interconnecting
microvilli in rafts; cells (green spheres, not to scale) are at the vertices
of a geodesic icosahedron. Blue arrows indicate flagella forces, red segments
represent microvilli, blue dots the microvilli tips, and blue lines the
collar-collar interface. (b) Connectivity of the whole raft. (c) Connections
between two cells. Effect of curvature $1/R$ on (d) geometry of raft and (e)
swimming speed $U$ (blue) and flow rate $Q$ (red) passing through $S_{f}$ at
constant. (f,g) Cross section of the disturbance flows
$\tilde{\mathbf{u}}_{d}$ and $\mathbf{u}_{f}$ around the a raft ($\Phi=0.31$,
$R=8.61$) in the reciprocal problems for calculating $U$ and $Q$.

As a first step toward understanding principles that govern the behavior of
such a novel organism as C. flexa, we analyze two models for these shape-
shifting structures. First, the fluid mechanics are studied by representating
the cell raft as a collection of spheres distributed on a hemispherical
surface, with nearby point forces to represent the action of flagella. Such a
model has been used to describe the motility of small sheet-like multicellular
assembles such as the alga Gonium GoniumPRE . The motility and filtering flow
through these rafts as a function of cell spacing and curvature explain the
observed properties of C. flexa. Second, abstracting the complex elastic
interactions between cells to the simplest connectivity, we show that a model
based on linear elasticity at the microscale produces bistability on the
colony scale.

_Fluid mechanics of feeding and swimming_ – The cells in a C. flexa raft are
ellipsoidal, with major and minor axes $a\sim 4\,\mu$m and $b\sim 3\,\mu$m,
with a single flagellum of length $2L\sim 24.5\,\mu$m and radius $r\sim
0.5\,\mu$m beating with amplitude $d\sim 2.3\,\mu$m and frequency
$f\sim\\!44\,$Hz SM , sending bending waves away from the body. A cell swims
with flagellum and collar rearward; the body and flagellum comprise a “pusher”
force dipole. From resistive force theory Laugabook  we estimate the flagellar
propulsive force to be $F\sim
2L\left(\zeta_{\perp}-\zeta_{\parallel}\right)(1-\beta)f\lambda\sim 8\,$pN,
where $\beta$ is a function of the wave geometry, $\lambda\sim 12\,\mu$m is
the wavelengths along the direction of the bending wave SM , $\zeta_{\perp}$
and $\zeta_{\parallel}$ are transverse and longitudinal drag coefficients,
$\zeta_{\perp}\sim 2\zeta_{\parallel}\sim 4\pi\mu/\ln(2L/r)$, with $\mu$ the
fluid viscosity. These features motivate a computational model in which $N$
identical cells in a raft have a spherical body of radius $a$ and a point
force $F\hat{\bf n}_{i}$ acting on the fluid a distance $L$ from the sphere
center, oriented along the vector $\hat{\bf n}_{i}$ that represents the collar
axis [Fig. 2(a,b)]. An idealization of a curved raft involves placing those
spheres on a connected subset of the vertices of a geodesic icosahedron (one
whose vertices lie on a spherical surface) of radius $\rho\gg a$; the area
fraction $\Phi$ of the sheet occupied by cells scales as $\Phi\sim
N(a/\rho)^{2}$. The pentagonal neighborhoods within the geodesic icosahedron
serve as topological defects that allow for smooth large-scale surface
curvature seung1988defects . Importantly, confocal imaging of C. flexa
colonies shows that a significant fraction ($\sim\\!\\!0.25$) of the cellular
neighborhoods defined by the microvilli connections are pentagonal SM , and
earlier work on C. perplexa leadbeater1983life  also found non-hexagonal
packing. We use the geodesic icosahedron $\\{3,5+\\}_{(3,0)}$ in standard
notation Wenninger , with $92$ total vertices, and take patches with $N=58$
for computational tractability. The vectors $\hat{\bf n}_{i}$ point towards
(away from) the icosahedron center in the flag-in (flag-out) forms (Fig.
2(a)). A deformation of the sheet to a new radius $\rho^{\prime}$, at fixed
$\Phi$, requires the new polar angle $\theta_{i}^{\prime}$ of a cell with
respect to the central axis of the sheet be related to its original angle
$\theta_{i}$ via
${\rho^{\prime}}^{2}\left(1-\cos{\theta_{i}}^{\prime}\right)=\rho^{2}\left(1-\cos\theta_{i}\right)$.
We define the scaled force offset length $\ell=L/a\sim 3$ and sheet radius
$R=\rho^{\prime}/a\gtrsim 6$, and take $R>0$ in the flag-in state.

Images of many colonies of C. flexa SM  show that the packing fraction in the
flag-out state $\Phi_{\rm out}=0.47\pm 0.06$, considerably less than both the
maximum packing fraction $\Phi_{\rm max}=\pi\sqrt{3}/6\simeq 0.907$ for a
hexagonal array of spheres in a plane, and the estimated maximum packing
fraction $\tilde{\Phi}_{\rm max}\simeq 0.83$ for circles on a sphere
circlepacking . The packing fraction in the flag-in state is $\Phi_{\rm
in}=0.34\pm 0.03$, and we use the extremes $\Phi_{\rm out}=0.53$ and
$\Phi_{\rm in}=0.31$ as representative values to explore the consequences of
the differences between the two forms.

Consider first an isolated force-free spherical cell at the origin moving at
velocity $U_{s}\hat{\mathbf{e}}_{x}$ with point force $-{\bf
F}=-F\hat{\mathbf{e}}_{x}$ at $-L\hat{\mathbf{e}}_{x}$ acting on the fluid and
its reaction force ${\bf F}$ acting on the cell. The cell experiences Stokes
drag $-\zeta_{s}U_{s}\hat{\mathbf{e}}_{x}$, where $\zeta_{s}=6\pi\mu a$, and a
disturbance drag $-D\hat{\mathbf{e}}_{x}$ arising from the disturbance flow
created by the point force. By the reciprocal theorem RT , the disturbance
drag is ${D}={\bf
F}\cdot\tilde{\mathbf{u}}_{d}(-L\hat{\mathbf{e}}_{x})/\hat{U}$, where
$\tilde{\mathbf{u}}_{d}({\bf r})$ is the disturbance flow created when the
cell is dragged along $\hat{\mathbf{e}}_{x}$ with unit speed $\hat{U}$. Force
balance then yields the single-cell swimming speed
$U_{s}\equiv(F/\zeta_{s})[1-3/(2\ell)+1/(2\ell^{3})]$. Thus, the closer the
point force is to the cell (i.e., the smaller is $\ell$), the more drag the
cell experiences and the slower is ${U}_{s}$. Setting $\ell=3$ yields
$U_{s}\sim 118\,\mu$m/s, consistent with observations.

This intuitive picture extends to a raft of cells. As the raft moves at
velocity $U\hat{\mathbf{e}}_{x}$, it experiences a Stokes drag $-\zeta
U\hat{\mathbf{e}}_{x}$. The disturbance flow created by the point forces
$F\hat{\bf n}_{i}$ acting at $\mathbf{r}_{i}+L\hat{\mathbf{n}_{i}}$, produces
a disturbance drag $D=F\Sigma_{i}\hat{\bf
n}_{i}\cdot\mathbf{u}_{d}(\mathbf{r}_{i}+L\hat{\mathbf{n}_{i}})$, where
${\mathbf{u}}_{d}$ is the (dimensionless) disturbance flow from the raft when
it is dragged along $\hat{\mathbf{e}}_{x}$ with unit speed. Force balance then
yields

$U=-\frac{F}{\zeta}\sum_{i=1}^{N}\hat{\bf
n}_{i}\cdot\left[\hat{\mathbf{e}}_{x}-{\mathbf{u}}_{d}\left(\mathbf{r}_{i}+L\hat{\bf
n}_{i}\right)\right],$ (1)

where $-F\Sigma_{i}\hat{\bf n}_{i}$ is the sum of reaction forces propelling
the raft along $\hat{\mathbf{e}}_{x}$, and where ${\mathbf{u}}_{d}$ has been
rendered dimensionless by the unit speed. In practice, we compute
${\mathbf{u}}_{d}$ and $\zeta$ using a Boundary Element Method GoniumPRE .
Because of the curved geometry, point forces are closer to neighboring cells
in the flag-in state than in the flag-out state. Thus, as in Fig. 2(d), for a
geometry with a given $|R|$, the flag-in state has a larger disturbance drag
than the flag-out state, and a smaller speed $U$.

The difference in swimming speed between the two states can also be explained
in terms of ${\mathbf{u}}_{d}(\mathbf{r}_{i}+L\hat{\bf n}_{i})$ in (1). Figure
2(f) shows that ${\mathbf{u}}_{d}$ inside the raft is close to
$\hat{\mathbf{e}}_{x}$ because of the curved geometry and screening effects.
Hence, $U$ is small when the point forces are inside. Meanwhile,
${\mathbf{u}}_{d}$ outside decays with the distance from the raft, so $U$ is
large when $\mathbf{r}_{i}+L\hat{\bf n}_{i}$ is outside.

Previous work on filter-feeding in choanoflagellates focused first on the far-
field limit based on a stresslet description Roper , but later work showed
near-field effects can significantly affect capture rates Kirkegaard . To
estimate the filter-feeding flux $Q$ passing through a colony of C. flexa, we
measure, in the body frame, the flux passing through the surface $S_{f}$
projected a distance of $1.2a$ from the cell center along $\hat{\bf n}$, as in
Fig. 2(g). By the reciprocal theorem, $Q$ can be written in terms of the
disturbance flow ${\mathbf{u}}_{f}$ around a stationary raft and the
hydrodynamic forces ${\mathbf{F}}_{f}$ on the raft when the surface $S_{f}$
applies a unit normal pressure $\hat{\mathbf{p}}$ on the fluid,

$Q=F\sum_{i}\hat{\mathbf{n}}_{i}\cdot{\mathbf{u}}_{f}(\mathbf{r}_{i}+L\hat{\bf
n}_{i})+{U}\hat{\mathbf{e}}_{x}\cdot\left({\mathbf{F}}_{f}-\int_{S_{f}}\hat{\mathbf{p}}\,dA\right),$
(2)

where ${\mathbf{u}}_{f}$ and ${\mathbf{F}}_{f}$ acquire the units of
velocity/pressure and area, respectively, by scaling with
$|\hat{\mathbf{p}}|$. Numerical results in Fig. 2(d) show that the flux due to
point forces $\sum_{i}F\hat{\mathbf{n}}_{i}\cdot{\mathbf{u}}_{f}/\hat{p}$
strongly dominates $Q$. Therefore, the difference in $Q$ between the two
states can be explained by ${\mathbf{u}}_{f}(\mathbf{r}_{i}+L\hat{\bf n}_{i})$
(Fig. 2(g)). To maintain incompressibility under pressure $\hat{\mathbf{p}}$,
the disturbance flow ${\mathbf{u}}_{f}$ is much stronger inside the raft than
outside. Hence, point forces placed inside the raft pump more flow through the
raft than when placed outside.

Figure 2(d) shows the effect of changes in the raft curvature and packing
fraction. There is one $R$ that maximizes swimming speed in the flag-in state
and another one that maximizes feeding flux in the flag-out state. This arises
from a balance between the screening effect mentioned above and the alignment
of forcing. In the flag-out state, an initial decrease in curvature aligns the
forcing direction with the swimming direction, increasing swimming speed, but
a further reduction in curvature reduces the screening effect as cells are now
more spread out in the plane orthogonal to the swimming direction. A similar
argument applies to the flow rate maximum in the flag-in state. Comparing
these maxima, Fig. 2(d) shows that a spread-out colony results in more flux,
while a closely-packed colony results in faster motility. Thus, through the
interconversion between the two states, C. flexa takes advantage of the
hydrodynamics effect of the curved geometry for efficient filter-feeding and
swimming.

Figure 3: Inversion dynamics from numerical studies. (a)-(d) A colony,
initially at a hemispherical minimum with $(\phi_{0},\psi_{0})=(0.55,0.65)$,
inverts after a change to $(0.71,0.41)$, with $\ell_{0}=0.5$, $K_{\psi}=2$ and
$K_{\ell}=5$. Connections between collar vertices are shown in blue,
apicobasal axes as red arrows at cell body positions. (e) Minimum energy $E$
in the $\psi_{0}-\phi_{0}$ plane, where $C<0$ (flag-in) above the black line
$\psi_{0}=\phi_{0}$, and $C>0$ (flag-out) below. $(f)$ Evolution of $E$ vs $C$
as the colony relax towards a minimum energy state after instantaneous changes
in $(\phi_{0},\psi_{0})$ shown by the dotted and solid red and blue lines in
$(e)$.

_Mechanics of inversion_ – Detailed studies suggest that inversion requires an
active process within each cell, likely driven by contraction of an F-actin
ring at the apical pole through the action of myosin brunet2019light . Thus, a
full treatment would address the complex problem of elastic filaments
responding to the apical actomyosin system and adhering to each other. We
simplify this description by considering as in Fig. 2(e) that each cell $i$,
located at ${\bf r}_{i}$ and surrounded by $m_{i}$ neighbors, has $m_{i}$
rigid, straight filaments emanating from it. Two filaments from neighboring
cells $i$ and $j$ meet at vertex $\rho$ located at ${\bf r}_{\rho}$, with
$\phi_{i\rho}$ the angle between ${\bf r}_{\rho}-{\bf r}_{i}$ and the cell
normal vector $\hat{\bf n}_{i}$. Any two adjacent filaments emanating from
cell $i$, and which meet neighboring filaments at vertices $\rho$ and
$\sigma$, define a plane whose normal $\hat{\bf n}_{i\rho\sigma}$ points
toward the apicobasal axis $\hat{\mathbf{n}}_{i}$. That normal and its
counterpart $\hat{\bf n}_{j\sigma\rho}$ on cell $j$ determine the angle
$2\psi_{ij\rho\sigma}$ between the two planes.

As above, we use the geodesic icosahedron to define the cell positions and
thus determine the filament network connecting neighboring cells. The two sets
of angles $\\{\phi\\}$ and $\\{\psi\\}$ are used to define a Hookean elastic
energy that mimics the elasticity of the microvilli, allowing for preferred
intrinsic angles $\phi_{0}$ and $\psi_{0}$ that encode the effects of the
apical actomyosin system on the microvilli and the geometry of microvilli
adhesion. Allowing also for stretching away from a rest length $\ell_{0}$, the
energy is

$\displaystyle
E=\frac{1}{2}k_{\phi}\sum_{i,\rho}\delta\phi_{i\rho}^{2}+\frac{1}{2}k_{\psi}\\!\\!\sum_{i,j,\rho,\sigma}\\!\\!\delta\psi_{ij\rho\sigma}^{2}+\frac{1}{2}k_{\ell}\sum_{i,\rho}\delta\ell_{i\rho}^{2},$
(3)

where $\delta\phi_{i\rho}=\phi_{i\rho}-\phi_{0}$,
$\delta\psi_{ij\rho\sigma}=\psi_{ij\rho\sigma}-\psi_{0}$, and
$\ell_{i\rho}=|{\bf r}_{i}-{\bf r}_{\rho}|-\ell_{0}$. The moduli $k_{\phi}$,
$k_{\psi}$ and $k_{\ell}$ and quantities $\phi_{0},\psi_{0}$ and $\ell_{0}$
are assumed constant for all cells.

The energy (3) is intimately tied to the lattice geometry of the raft. If the
cells are arranged in a hexagonal lattice ($m_{i}=6$) the system of filaments
can achieve $E=0$ by setting all cell-collar angles to $\phi_{0}$, all collar-
collar interface angles to $\psi_{0}$, and $\phi_{0}=\psi_{0}$. This
corresponds to a flat sheet. Increasing $\phi_{0}=\psi_{0}$ leads to uniform,
isotropic sheet expansion. In a non-planar raft, curvature is introduced
through topological defects ($m_{i}\neq 6$), such as pentagons, and mismatch
between the local values of $\psi$ and $\phi$. While pentagonal defects are
known to cause out-of-plane buckling in crystal lattices seung1988defects ,
they do not by themselves select a particular sign of the induced curvature.
Thus, there is inherent bistability in the cellular raft that can be biased by
changes in the geometry of the out-of-plane filaments, somewhat akin to the
role of “apical constriction” in the shapes of epithelia Hannezo .

For the case of two cells lying in a plane, each with two filaments, and with
one vertex between them, if $\phi=\phi_{0}$, $\psi=\psi_{0}$, and
$r=\ell_{0}$, then the filament tips lie on a circle of radius
$R_{0}=1/C_{0}$, where $C_{0}=\sin(\psi_{0}-\phi_{0})/\ell_{0}\sin\phi_{0}$.
While, in general, the equilibrium state of a curved raft will not have
$\phi_{i\rho}=\phi_{0}$, $\psi_{ij\rho\sigma}=\psi_{0}$ and
$r_{i\rho}=\ell_{0}$ everywhere, we may nevertheless use this relationship to
define a proxy for the average curvature of the raft. Recognizing that in
numerical studies stretching effects are small, we ignore variations in
$r_{i\rho}$ and define
$C=\sin(\langle\psi\rangle-\langle\phi\rangle)/\ell_{0}\sin(\langle\phi\rangle)$,
where $\langle\cdot\rangle$ is an average over cells and vertices. The colony
is in the flag-in (flag-out) state when $C>0$ ($C<0$).

The simplest model of raft dynamics localizes the viscous drag to the
individual cell and vertex positions ${\bf r}_{\gamma}$ according to a
gradient flow $\zeta\partial_{t}{\bf r}_{\gamma}=-\partial
E/\partial\mathbf{r}_{\gamma}$ driven by the force derived from (3). We solve
this dynamics numerically with forward integration. Since the
$\hat{\mathbf{n}}_{i}$ are constrained to have unit length, they are
normalized after each step in the direction of the negative gradient, making
the dynamical algorithm follow a projected gradient descent eicke1992iteration
. Via a rescaling of time we may set one of the elastic constants to unity
(say, $k_{\phi}$) and need only consider the ratios
$K_{\psi}=k_{\psi}/k_{\phi}$ and $K_{\ell}=k_{\ell}/k_{\phi}$.

Interconversion between the flag-in and flag-out states is shown in Figs.
3(a-d) following an abrupt change in the preferred angle pair
$(\phi_{0},\psi_{0})$ that crossing the line of equality $\psi_{0}=\phi_{0}$
that divides the states (Fig. 3(e)), as during raction/relaxation of the
F-actin ring in response to a stimulus. The intermediate shapes exhibit a ring
of inflection points similar to those seen in experiments on C. flexa and also
in the inversion the algae Pleodorina hohn2016distinct  and larger species
viamontes1977cell ; PLOS . Tracking the energy as each of the two equilibria
is achieved, the picture that emerges in Fig. 3(f) is evolution on a double-
well potential energy landscape as a biasing field is switched in sign.

We have shown that simple models can explain the swimming, feeding, and
inversion of the recently discovered multicellular choanoflagellate C. flexa
brunet2019light . These results suggest further exploration on a possible
continuum description of the sheets, fluid-structure interactions during
locomotion, dynamics of photokinesis, and developmental processes of these
remarkable organisms.

###### Acknowledgements.

We gratefully acknowledge Gabriela Canales and Tanner Fadero for high speed
imaging assistance and Kyriacos Leptos for comments and suggestions. This work
was supported in part by a Research Fellowship from Peterhouse, Cambridge
(LF), a Churchill Scholarship (AK), JSPS Kakenhi (TI), The John Templeton
Foundation and Wellcome Trust Investigator Grant 207510/Z/17/Z (REG). For the
purpose of open access, the authors have applied a CC BY public copyright
license to any Author Accepted Manuscript version arising from this
submission.

## References

  * (1) L. Solnica-Krezel and D.S. Sepich, Gastrulation: Making and Shaping Germ Layers, Annu Rev. Cell Dev. Biology 28, 687–717 (2012).
  * (2) A. Goriely, D.E. Moulton and R. Vandiver, Elastic cavitation, tube hollowing, and differential growth in plants and biological tissues, EPL 91, 18001 (2010).
  * (3) S. Höhn and A. Hallmann, Distinct shape-shifting regimes of bowl-shaped cell sheets–embryonic inversion in the multicellular green alga Pleodorina, BMC Dev. Biol. 16, 35 (2016).
  * (4) S. Höhn, A.R. Honerkamp-Smith, P.A. Haas, P. Khuc Trong, and R.E. Goldstein, Dynamics of a Volvox embryo turning itself inside out, Phys. Rev. Lett. 114, 178101 (2015).
  * (5) T. Brunet and N. King, The Origin of Animal Multicellularity and Cell Differentiation, Dev. Cell 43, 124–140 (2017).
  * (6) S.R. Fairclough, M.J. Dayel, and N. King, Multicellular development in a choanoflagellate, Curr. Biology 20, R875–R876 (2010).
  * (7) B.T. Larson, T. Ruiz-Herrero, S. Lee, S. Kumar, L. Mahadevan and N. King, Biophysical principles of choanoflagellate self-organization, Proc. Natl. Acad. Sci. USA 117, 1303–1311 (2020).
  * (8) T. Brunet, B.T. Larson, T.A. Linden, M.J.A. Vermeij, K. McDonald, and N. King, Light-regulated collective contractility in a multicellular choanoflagellate, Science 366, 326–33 (2019).
  * (9) As noted previously brunet2019light , the behavior of C. flexa is similar to the species C. perplexa studied earlier: B.S.C. Leadbeater, Life-history and ultrastructure of a new marine species of Proterospongia (Choanoflagellida), J. Mar. Biol. Ass. U.K. 63, 135–160 (1983).
  * (10) H. de Maleprade, F. Moisy, T. Ishikawa, and R.E. Goldstein, Motility and phototaxis in Gonium, the simplest differentiated colonial alga Phys. Rev. E 101, 022416 (2020).
  * (11) See Supplemental Material at http://link.aps.org/supplemental/xxx for further experimental results, which includes Ref. [12].
  * (12) J. Schindelin, I. Arganda-Carreras, E. Frise et al. Fiji: an open-source platform for biological-image analysis. Nat Methods 9, 676–682 (2012).
  * (13) E. Lauga, The Fluid Dynamics of Cell Motility (Cambridge University Press, Cambridge, UK, 2020).
  * (14) H.S. Seung and D.R. Nelson, Defects in flexible membranes with crystalline order, Phys. Rev. A 38, 1005–1018 (1988).
  * (15) M. Wenninger, Spherical Models (Cambridge University Press, Cambridge, UK, 1979).
  * (16) B.W. Clare and D.L. Kepert, The Optimal Packing of Circles on a Sphere, J. Math. Chem. 6, 325–349 (1991).
  * (17) S. Kim and S. J. Karrila, Microhydrodynamics: Principles and Selected Applications (Butterworth-Heinemann, Boston, US, 1991)
  * (18) Roper, M. and Dayel, M.J. and Pepper, R.E. and Koehl, M.A.R., Cooperatively Generated Stresslet Flows Supply Fresh Fluid to Multicellular Choanoflagellate Colonies, Phys. Rev. Lett. 110, 228104 (2013).
  * (19) J.B. Kirkegaard and R.E. Goldstein, Filter-feeding, near-field flows, and the morphologies of colonial choanoflagellates, Phys. Rev. E 94, 052401 (2016).
  * (20) E. Hannezo, J. Prost, J.-F. Joanny, Theory of epithelial sheet morphology in three dimensions, Proc. Natl. Acad. Sci. USA 111, 27–32 (2014).
  * (21) B. Eicke, Iteration methods for convexly constrained ill-posed problems in Hilbert space Num. Func. Anal. Opt. 13, 413–429 (1992).
  * (22) G.I. Viamontes and D.L. Kirk, Cell shape changes and the mechanism of inversion in Volvox, J. Cell Biol. 75, 719–730 (1977).
  * (23) P.A. Haas, S. Höhn, A.R. Honerkamp-Smith, J.B. Kirkegaard, and R.E. Goldstein, The Noisy Basis of Morphogenesis: Mechanisms and Mechanics of Cell Sheet Folding Inferred from Developmental Variability, PLOS Biol. 16, e2005536 (2018).

##  I Supplemental Material

This file contains additional experimental results on flagellar dynamics and
geometry of C. flexa.

I. VIDEO IMAGING AND ANALYSIS

A. Supplementary Videos

C. flexa sheets were imaged in FluoroDishes (World Precision Instruments
FD35-100) by differential interference contrast (DIC) microscopy using a
$40\times$ (water immersion, C-Apochromat, 1.1 NA) Zeiss objective mounted on
a Zeiss Observer Z.1 with a pco.dimax cs1 camera.

Flagellar characteristics reported in Table 1 were obtained as follows. Beat
frequencies were determined by averaging over five cycles for each of twenty
randomly selected cells. All other measurements are averages over ten randomly
selected cells. The comparatively large wavelength for flag-out sheets may be
due in part to the fact that the flagellar waveform in that state is not
sinusoidal, and its wavelength is thus less well defined than in the flag-in
state.

Figure S1: Flagellar dynamics in C. flexa colonies. (a) Snapshot from Video 1, a high speed recording of a flag-in sheet used to determine flagellar characteristics reported in Table 1. Movie is set to play at $17\times$ slower than real time. (b) As in (a), but for the flag-out state in Video 2. Scale bars are $10\,\mu$m.  conformation | beat frequency $f$ | length $2L$ | amplitude $d$ | wavelength $\lambda$
---|---|---|---|---
flag-in |  $45\pm 4$ Hz |  $26\pm 5$ $\mu$m |  $2.4\pm 0.6$ $\mu$m |  $9\pm 3$ $\mu$m
flag-out |  $43\pm 8$ Hz |  $23\pm 4$ $\mu$m |  $2.2\pm 0.7$ $\mu$m |  $15\pm 3$ $\mu$m
Table 1: Measurements of flagellar characteristics for the flag-in and flag-
out sheets in Videos S1 and S2. Uncertainties reported are standard
deviations.

B. Estimating the propulsive force from the flagella

In the fluid mechanics model of the C. Flexa raft, we approximate the flagella
beating as an effective propulsive point force $\mathbf{F}$ acting in the
direction $\hat{\mathbf{n}}$. The magnitude of this force $F$ can be
approximated using the resistive force theory Laugabook1  as

$F=2L\left(\zeta_{\perp}-\zeta_{\parallel}\right)(1-\beta)f\lambda$ (S1)

where $2L$ is the flagella length, $f$ the beat frequency and $\lambda$ the
projected wavelength in the direction of the traveling sinusoidal wave (i.e.
$\hat{\mathbf{n}}$), the values of which are listed in Table 1. Meanwhile,

$\zeta_{\perp}=\frac{4\pi\mu}{\ln(2L/r)}\quad\mbox{and}\quad\zeta_{\parallel}=\frac{2\pi\mu}{\ln(2L/r)}$
(S2)

are the transverse and longitudinal drag coefficients of a cylindrical
filament of radius $r$, approximated using the resistive force theory, and
$\beta$ is a coefficient that depends on the flagella waveform. Although the
flagella waveform is not necessarily sinusoidal, in the absence of better
measurements, the value of $\beta$ is approximated, assuming the flagella
takes a sinusoidal waveform $f(x)$ with wavelength $\lambda$ and amplitude $d$
(Table 1), as

$\beta=\frac{\int_{0}^{\lambda}\left(1+f^{\prime}(x)\right)^{-1/2}dx}{\int_{0}^{\lambda}\left(1+f^{\prime}(x)\right)^{1/2}dx},$
(S3)

which can be found by numerically. In the limiting of $2\pi d/\lambda\ll 1$,
$\beta\approx 2\pi^{2}(d/\lambda)^{2}$.

II. CONFOCAL IMAGING

Sheets in Figs. S2 and S3 were fixed and stained with FM1-43FX or with Alexa
488-phalloidin as in brunet2019light1 . Sheets were imaged on Zeiss LSM 880
with AiryScan using a 63x, 1.4 NA C Apo oil immersion objective (Zeiss).
Z-projections were generated with Fiji schindelin . Packing fraction was
estimated by projecting cell bodies located within the same plane in a locally
flat portion of the sheet and by manually outlining the border of the colonies
(red dotted line in Fig. S2).

Figure S2: Packing fraction in flagella-in and flagella-out C. flexa colonies.
(a) 3D stacks (left column) of sheets stained with the fluorescent membrane
marker FM 1-43 FX were imaged by confocal microscopy to determine sheet
morphology. Packing fraction was computed by doing a Z-projection of a locally
flat portions of individual sheets (middle column) and generating a binarized
image in which the area occupied by individual cells appears black (right
column). Packing fraction is the ratio of the area occupied by cells to the
total colony area within that plane (area within the red dotted line). (b)
Boxplot depicting packing fraction values for 6 sheets with flagella out and 6
sheets with flagella-in. p=0.2% by the Mann-Whitney test.

Packing fraction was then computed as the ratio of the area occupied by cells
(black area in binarized image in Fig. S2) to the total area occupied by the
colony (area within the red dotted line in Fig. S2). Polygonal collar borders
in Fig. S3 were manually outlined and colored with Adobe Illustrator 27.3.1
(2023). Hexagonal and pentagonal outlines were counted in 9 colonies and
counts are reported in Table 2.

Figure S3: Hexagonal and polygonal collar outlines within flagella-in and flagella-out C. flexa colonies. F-actin within colonies was stained with fluorescent phalloidin, which outlines the microvillous collars linking cells, as well as the actin cytoskeleton within the cell body. Left column: whole-colony Z-projections. Middle left: Z-projections of a few planes intersecting many collars, showing the polygonal outlines of collar contacts. Middle right: color-coded polygons showing a majority of hexagons and a minority of pentagons. Right: whole-colony Z-projection with the plane containing most collars outlined in purple. Table 2: Numbers of hexagons and pentagons in representative sections of $9$ colonies of C. flexa. colony # | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $8$ | $9$ | total
---|---|---|---|---|---|---|---|---|---|---
hexagons | 4 | 4 | 2 | 8 | 8 | 15 | 6 | 6 | 8 | 61
pentagons | 4 | 0 | 1 | 2 | 2 | 7 | 2 | 3 | 2 | 23

## References

  * (1) E. Lauga, The Fluid Dynamics of Cell Motility (Cambridge University Press, Cambridge, UK, 2020), §7.1.3-7.1.4.
  * (2) T. Brunet, B.T. Larson, T.A. Linden, M.J.A. Vermeij, K. McDonald, and N. King, Light-regulated collective contractility in a multicellular choanoflagellate, Science 366, 326–33 (2019).
  * (3) J. Schindelin, I. Arganda-Carreras, E. Frise et al. Fiji: an open-source platform for biological-image analysis. Nat Methods 9, 676–682 (2012).


# Some refinements for translation invariant
quadratic forms in dense sets

Eugen Keil  Mathematical Institute
University of Oxford
Andrew Wiles Building
Radcliffe Observatory Quarter
Woodstock Road
Oxford
OX2 6GG
United Kingdom<EMAIL_ADDRESS>
###### Abstract.

We improve the result of our previous paper on translation invariant quadratic
forms in two special cases. We reduce the density bound
$|\mathcal{A}|/N=O((\log\log N)^{-c})$ to $|\mathcal{A}|/N=O((\log N)^{-c})$
for most quadratic forms and handle almost diagonal equations in $s\geq 6$
variables instead of $s\geq 8$.

###### Key words and phrases:

dense set and quadratic form and translation invariant

###### 2010 Mathematics Subject Classification:

Primary 11B30; Secondary 11P55, 11D09

##  1\. Introduction

This paper complements the author’s work [3] on translation invariant
quadratic forms and should be read subsequent to it. In [3] we have proven the
following result.

######  Theorem 1.1.

Let $Q\in\mathbb{Z}^{s\times s}$ be symmetric with
$Q\cdot\boldsymbol{1}=\boldsymbol{0}$ and off-rank $r$. Assume that
$\mathbf{x}^{T}Q\mathbf{x}=0$ has a non-singular real solution and assume that
$s\geq 5+3r$ for $1\leq r\leq 4$ and $s\geq 10$ for $r\geq 5$. If there are
only trivial solutions, when the variables are restricted to
$\mathcal{A}\subset\\{1,2,\ldots,N\\}$, then $|\mathcal{A}|\leq
C_{Q}N(\log\log N)^{-c}$ for some $c,C_{Q}>0$ and $c$ independent of $Q$.

The condition $Q\cdot\boldsymbol{1}=\boldsymbol{0}$ is equivalent to
translation invariance of the equation $\mathbf{x}^{T}Q\mathbf{x}=0$ as
explained in [3]. The ‘off-rank’ of a matrix is defined in Section 2.

The aim of this note is to give improvements on this theorem in two natural
special cases. On the one hand, we can improve the bound $|\mathcal{A}|\leq
C_{Q}N(\log\log N)^{-c}$ to $|\mathcal{A}|\leq C_{Q}N(\log N)^{-c}$ in the
case $r\geq 5$ (which covers almost all quadratic forms) and on the other
hand, solve the problem with as few as $s\geq 6$ variables in the almost
diagonal situation $r=1$. In the end, we want to discuss how the two
approaches might be useful to improve on Theorem 1.1 for all quadratic forms.

The next section is devoted to explain the results, which are given in Theorem
2.2 and Theorem 2.4.

We introduce some standard notation. We use the Vinogradov notation $\ll$ and
$O$-notation throughout the paper and indicate dependencies on parameters by
subscripts, like in $O_{P,\epsilon}(N)$. The asymptotic parameter
$N\in\mathbb{N}$ should be thought of as large and we work most of the time on
the interval $[1,N]=\\{1,2,\ldots,N\\}$. We write
$\mathbb{T}=\mathbb{R}/\mathbb{Z}$ for the ‘circle’ and identify it with
$[0,1]$ whenever convenient. As usual we abbreviate $e(x)=\exp(2\pi ix)$.

For a set $\mathcal{A}\subset[1,N]$ with the indicator function
$1_{\mathcal{A}}$ and density $\delta=|\mathcal{A}|/N$ we define the ‘balanced
function’ by

(1.1) $\displaystyle f(n)=1_{\mathcal{A}}(n)-\delta.$

We are using bold face to denote vectors and the notation $\mathbf{x}\leq N$
means that each coordinate is bounded by $N$.

Acknowledgements:
This work has a non-empty intersection with a chapter from the author’s Ph.D.
thesis. He would like to thank Trevor Wooley for his encouragement and helpful
discussions. The Ph.D.-studies of the author were partially supported by the
EPSRC. This work was finished at the University of Oxford with the support by
EPSRC, grant EP/J009458/1.

##  2\. Statement of Results

First we summarize some of the definitions from [3] that are needed for the
proofs below.

A quadratic form $Q(\mathbf{x})=\mathbf{x}^{T}Q\mathbf{x}$ is translation
invariant if $Q(\mathbf{x}+\boldsymbol{1})=Q(\mathbf{x})$ for all
$\mathbf{x}\in\mathbb{Z}^{s}$. For a corresponding symmetric matrix
$Q\in\mathbb{Z}^{s\times s}$, this translates into the condition
$Q\cdot\boldsymbol{1}=\boldsymbol{0}$.

######  Definition 2.1 (Off-diagonal rank).

For a symmetric matrix $Q\in\mathbb{R}^{s\times s}$ we consider the set of
matrices $M$ such that $M=P^{T}QP$ for a permutation matrix $P$. Such a matrix
can be (non-uniquely) written as

$\displaystyle M=\begin{pmatrix}A&B\\\ B^{T}&C\end{pmatrix}$

for some matrices $A,B$ and $C$. The _off-rank_ $r$ of $Q$ is defined as

$\displaystyle r=\max\mbox{rank}(B),$

where the maximum is taken over all choices of $P$ and decompositions of $M$.
In other words, $r$ is the maximal rank of a submatrix in $Q$, that doesn’t
contain any diagonal elements.

Now we can state the first theorem.

######  Theorem 2.2.

Let $Q\in\mathbb{Z}^{s\times s}$ be symmetric with
$Q\cdot\boldsymbol{1}=\boldsymbol{0}$ and off-rank $r\geq 5$ (this implies
that $s\geq 10$). Assume that it has a non-singular real solution to
$\mathbf{x}^{T}Q\mathbf{x}=0$, but only trivial solutions when the variables
are restricted to $\mathcal{A}\subset[1,N]$. Then $|\mathcal{A}|\leq
C_{Q}N(\log N)^{-c}$ for some $c,C_{Q}>0$.

The main idea for the proof can be summarized as follows. The usual density
increment procedure depends on the fact that $e(\alpha p(x))$ for a polynomial
$p$ is constant on long arithmetic progressions. If $\alpha$ is close to a
rational number with small denominator (on the ‘major arcs’), the length of
those progressions is comparable to $N$ and we get a much better bound. It
turns out that we can control the contribution of the ‘minor arcs’ (complement
of the major arcs) in the proof of Theorem 2.2 in a uniform way independent of
the structure of $\mathcal{A}$. This is done in Section 3. In Section 4 we
perform the density increment with $\alpha$ in the major arcs.

To motivate the second part of the paper, consider quadratic forms defined by
the following matrices.

######  Example 2.3.

Consider the matrices

$\displaystyle\begin{pmatrix}3&1&0&0&0\\\ 1&7&-2&-3&2\\\ 0&-2&5&0&0\\\
0&-3&0&1&0\\\ 0&2&0&0&-4\end{pmatrix}\qquad\mbox{ and
}\qquad\begin{pmatrix}0&2&2&2&1\\\ 2&-1&4&4&2\\\ 2&4&6&4&2\\\ 2&4&4&-5&2\\\
1&2&2&2&3\end{pmatrix}.$

The first example is ‘almost diagonal’ in the sense that off-diagonal terms
are all concentrated in the second column/row. In what sense is the second
example almost diagonal? We can write the second matrix as a sum of a diagonal
matrix $D$ and a rank one perturbation $R$, a matrix of the form
$R=\mathbf{v}\cdot\mathbf{v}^{T}$ for the vector $\mathbf{v}=(1,2,2,2,1)^{T}$.
While the first matrix is an example of a small ‘local’ perturbation, the
second is an instance of a diagonal matrix with a small ‘global’ perturbation.

One can easily check that the two matrices have off-rank one. As we will see
in Section 5, it turns out that those two cases exhaust the possible ways that
a symmetric matrix can have off-rank one. Armed with this classification, we
can show that Theorem 1.1 holds for $s\geq 6$ (instead of $s\geq 5+3r=8$) if
the matrix underlying our quadratic form has off-rank one. The result is given
in the following theorem.

######  Theorem 2.4.

Let $Q\in\mathbb{Z}^{s\times s}$ be symmetric with
$Q\cdot\boldsymbol{1}=\boldsymbol{0}$ and off-rank $r=1$. Assume that $s\geq
6$, and that there is a non-singular real solution to
$\mathbf{x}^{T}Q\mathbf{x}=0$. If there are only trivial solutions when the
variables are restricted to $\mathcal{A}\subset[1,N]$, then
$|\mathcal{A}|\ll_{Q}N(\log\log N)^{-1/15}$.

Let us analyse the two cases in more detail. What does it mean that in the
first example all off-diagonal terms are concentrated in the row/column of one
variable? Let us assume for simplicity that the variable is $x_{s}$. We will
see that by completing squares, we can remove all off-diagonal terms and
replace each variable $x_{i}$ by $x_{i}-x_{s}$ (translation invariance!),
reducing the problem to a diagonal form in five variables, which is well
within reach of the classical approach. We deal with this case in Section 6.

In Section 7 we reduce the second case to the following diagonal system
considered by Smith [5] and the author [2].

(2.1) $\begin{split}d_{1}x_{1}^{2}+d_{2}x_{2}^{2}+\ldots+d_{s}x_{s}^{2}&=0,\\\
d_{1}x_{1}+d_{2}x_{2}+\ldots+d_{s}x_{s}&=0.\end{split}$

The system (2.1) can be handled under the following conditions.

######  Conditions for system (2.1).

$(i)$ $d_{1}+d_{2}+\ldots+d_{s}=0$,
$(ii)$ $s\geq 7$ and $d_{i}\neq 0$ for all $1\leq i\leq s$,
$(iii)$ there are at least two positive and at least two negative coefficients
$d_{i}$.

Condition $(i)$ encodes the translation invariance of the system. The bound
$s\geq 7$ in condition $(ii)$ will correspond to the condition $s\geq 6$ in
Theorem 2.4. We assume that $d_{i}\neq 0$ since we can easily locate non-
trivial solutions otherwise. Condition $(iii)$ is needed for the existence of
a non-singular real solution (see [2] for details).

In Section 7 we use the following theorem from [2] to deduce Theorem 2.4 in
the second case.

######  Theorem 2.5.

Assume that the conditions above hold and system (2.1) has only trivial
solutions for $x_{i}\in\mathcal{A}\subset[1,N]$. Then $|\mathcal{A}|\leq
CN(\log\log N)^{-1/15}$ for some constant $C$, which depends only on the
coefficients $d_{i}$ of the system.

##  3\. Controlling the minor arcs

As in [3], the main ingredient in the proof is a bilinear sum estimate for our
exponential sum

$\displaystyle S_{g}(\alpha)=\sum_{\mathbf{x}\leq N}g(\mathbf{x})e(\alpha
Q(\mathbf{x})).$

Write $\alpha=a/q+\beta$ for a diophantine approximation with $q\leq
N,|\beta|\leq(qN)^{-1}$ and

(3.1) $\displaystyle K(\alpha)=\Big{(}N\log
q+\min\Big{\\{}\frac{N^{2}}{q},\frac{|\log(|\beta|N^{2})|+1}{|\beta|q}\Big{\\}}\Big{)}^{1/2}.$

Then we can bound $S_{g}(\alpha)$ as in [3] for $|g|\leq 1$ by

$\displaystyle|S_{g}(\alpha)|\ll N^{s-r}K(\alpha)^{r}.$

As in [3, Theorem 3], we can refine the estimate if we are given
$L^{1}$-bounds $\sum_{n\leq N}|g(n)|=O(\delta N)$. As in [3], we get the bound

(3.2) $\displaystyle|S_{g}(\alpha)|\ll\delta^{s-10}N^{s-5}K(\alpha)^{5}.$

This pointwise bound allows us to deduce a sharp $L^{p}$-estimate for
$S_{g}(\alpha)$, which is necessary for the density increment strategy to
work.

If we look at the proof of this $L^{p}$ bound in [3], we see that the
condition we need for the exponent of $K(\alpha)$ is that it is bigger than
four. This means that we can pull out a small power of $K(\alpha)$, and
improve the estimate on the minor arcs, where $K(\alpha)$ is small.

For an absolute constant $D>1$ to be chosen later, define the major arcs for
$q\leq D^{4}\delta^{-40}$ and $(a;q)=1$ by

(3.3)
$\displaystyle{\mathfrak{M}}(q,a)=\\{\alpha\in\mathbb{T}:\|\alpha-a/q\|\leq
D^{8}\delta^{-80}N^{-2}\\}$

(disjoint for $N\geq D^{8}\delta^{-80}$) and set ${\mathfrak{M}}$ to be the
union of all those sets. The minor arcs
${\mathfrak{m}}=\mathbb{T}\backslash{\mathfrak{M}}$ is the complement. The key
point to notice is that the constants in the definition of ${\mathfrak{M}}$
only depend on $\delta$ and an absolute constant $D$. The precise numbers are
less important.

By Dirichlet’s approximation theorem, we have that every $\alpha\in\mathbb{T}$
is contained in at least one ball defined by $\|\alpha-a/q\|\leq(qN)^{-1}$ for
some $q\leq N$. This means that we have two types of minor arcs. Those
Dirichlet neighbourhoods with $q>D^{4}\delta^{-40}$ and those with $q\leq
D^{4}\delta^{-40}$ but $D^{8}\delta^{-80}N^{-2}<\|\alpha-a/q\|\leq(qN)^{-1}$.

If $q>D^{4}\delta^{-40}$, we have

$\displaystyle K(\alpha)\ll(N\log q)^{1/2}+\frac{N}{q^{1/2}}\ll
D^{-2}\delta^{20}N.$

for $N\gg D^{8}\delta^{-80}$.

In the second case $q\leq D^{4}\delta^{-40}$, we split the two intervals
$D^{8}(N^{2}\delta^{80})^{-1}<|\beta|\leq(qN)^{-1}$ into dyadic bits
$|\beta|\in(2^{i}N^{-2},2^{i+1}N^{-2}]$ for $\log_{2}(\delta^{-80}D^{8})\leq
i\leq\log_{2}(Nq^{-1})$. For fixed $i$ we have

$\displaystyle K(\alpha)$ $\displaystyle\ll(N\log
q)^{1/2}+\Big{(}N^{2}\frac{i+1}{q2^{i}}\Big{)}^{1/2}$ $\displaystyle\ll(N\log
q)^{1/2}+N2^{-i/4}\ll D^{-2}\delta^{20}N$

for $N\gg D^{8}\delta^{-80}$.

Now that we have good pointwise bounds, we cite [3, Lemma 2], which gives us a
control of $K(\alpha)$ on average. For $p>4$ we have

$\displaystyle\int_{0}^{1}|K(\alpha)|^{p}\,d\alpha\ll N^{p-2}.$

We apply the two bounds with (3.2) and get

$\displaystyle\int_{{\mathfrak{m}}}|S_{g}(\alpha)|\,d\alpha\ll$
$\displaystyle(\delta^{s-10}N^{s-5})\sup_{\alpha\in{\mathfrak{m}}}|K(\alpha)|^{1/2}\int_{{\mathfrak{m}}}|K(\alpha)|^{9/2}\,d\alpha$
$\displaystyle\ll$
$\displaystyle(\delta^{s-10}N^{s-5})(D^{-1}\delta^{10}N^{1/2})N^{5/2}\ll
D^{-1}\delta^{s}N^{s-2}.$

We can use this bound in the argument in [3, Section 3]. We know that

$\displaystyle\int_{0}^{1}|S_{g}(\alpha)|\,d\alpha\gg\delta^{s}N^{s-2}$

for $g(\mathbf{x})=1_{\mathcal{A}^{s}}(\mathbf{x})-\delta^{s}$ with
$\delta=|\mathcal{A}|/N$ as long as $N\gg_{Q}\delta^{-2}$. By choosing $D$
sufficiently small, this implies that

$\displaystyle\int_{{\mathfrak{M}}}|S_{g}(\alpha)|\,d\alpha\gg\delta^{s}N^{s-2}.$

By the same argument as in [3, Section 3] we get the bound

(3.4)
$\displaystyle\sup_{\alpha\in{\mathfrak{M}}}|S_{f_{i}}(\alpha)|\gg\delta^{s+80}N^{s}$

for a function $f_{i}$ that is a (tensor) product of functions $f,\delta$ and
$1_{\mathcal{A}}$ (see (1.1)).

##  4\. Density increment on major arcs

Consider the density increment argument from [3, Section 7]. Instead of a
general lower bound for the exponential sum, we have (3.4), where $\alpha$
lies in the major arcs ${\mathfrak{M}}$ defined in (3.3). This implies that
there is a $q\leq D^{4}\delta^{-40}$ and $1\leq a\leq q$ such that
$\|\alpha-a/q\|\leq D^{8}\delta^{-80}N^{-2}$. For the following argument, $D$
is a fixed constant and will be absorbed in the Vinogradov notation.

We decompose our interval $[1,N]$ into progressions modulo $q$. On each
progression we have $e(\alpha Q(\mathbf{x}))=e(\beta
Q(\mathbf{x}))e(aQ(\mathbf{x})/q)$, where $\alpha=a/q+\beta$. Since $\beta$ is
very small, we can estimate

$\displaystyle|e(\beta Q(\mathbf{x}))-e(\beta Q(\mathbf{y}))|\leq
2\pi|\beta||Q(\mathbf{x})-Q(\mathbf{y})|\ll(N^{2}\delta^{80})^{-1}|Q(\mathbf{x})-Q(\mathbf{y})|.$

The second factor is constant on progressions with difference $q$. If we write
$\mathbf{x}=q\mathbf{m}+\mathbf{h}$ and $\mathbf{y}=q\mathbf{n}+\mathbf{h}$,
we obtain

$\displaystyle|e(\alpha Q(\mathbf{x}))-e(\alpha Q(\mathbf{y}))|=|e(\beta
Q(\mathbf{x}))e(aQ(\mathbf{h})/q)-e(\beta Q(\mathbf{y}))e(aQ(\mathbf{h})/q)|$
$\displaystyle=$ $\displaystyle|e(\beta Q(\mathbf{x}))-e(\beta
Q(\mathbf{y}))|\ll(N^{2}\delta^{80})^{-1}|Q(\mathbf{x})-Q(\mathbf{y})|$
$\displaystyle\ll$
$\displaystyle(N^{2}\delta^{80})^{-1}q^{2}N\|\mathbf{m}-\mathbf{n}\|_{\infty}\ll
N^{-1}\delta^{-160}\|\mathbf{m}-\mathbf{n}\|_{\infty}.$

If we cut our $q$-progressions $P_{\mathbf{h}}$ into subprogressions
$P_{\mathbf{h},\mathbf{j}}$ of side length approximately $P=\eta\delta^{240}N$
for some small fixed $\eta$, our function $e(\alpha Q(\mathbf{x}))$ is
constant up to a small error term of size $O(\delta^{80})$. This implies

$\displaystyle\sum_{\mathbf{x}\leq N}f_{i}(\mathbf{x})e(\alpha
Q(\mathbf{x}))=\sum_{\mathbf{h}\leq
q}\sum_{\mathbf{j}\leq(q\eta\delta^{240})^{-1}}\sum_{\mathbf{x}\in
P_{\mathbf{h},\mathbf{j}}}f_{i}(\mathbf{x})e(\alpha Q(\mathbf{x}))$
$\displaystyle=$ $\displaystyle\sum_{\mathbf{h}\leq
q}\sum_{\mathbf{j}\leq(q\eta\delta^{240})^{-1}}\sum_{\mathbf{x}\in
P_{\mathbf{h},\mathbf{j}}}f_{i}(\mathbf{x})(c(\alpha,\mathbf{h},\mathbf{j})+\epsilon(\alpha,\mathbf{x})\delta^{80})$

for some functions $|c(\alpha,\mathbf{h},\mathbf{j})|\leq 1$ and
$|\epsilon(\alpha,\mathbf{x})|<1/2$. Taking absolute values and using 3.4 with
$\sum_{\mathbf{x}}|f_{i}(\mathbf{x})|\ll\delta^{s}N^{s}$ gives

$\displaystyle\sum_{\mathbf{h}\leq
q}\sum_{\mathbf{j}\leq(q\eta\delta^{240})^{-1}}\Big{|}\sum_{\mathbf{x}\in
P_{\mathbf{h},\mathbf{j}}}f_{i}(\mathbf{x})\Big{|}\gg\delta^{s+80}N^{s}.$

For at least one value of $\mathbf{h}$ and $\mathbf{j}$, we get

$\displaystyle\Big{|}\sum_{\mathbf{x}\in
P_{\mathbf{h},\mathbf{j}}}f_{i}(\mathbf{x})\Big{|}\gg\delta^{240s+s+80}N^{s}.$

Since $f_{i}$ is a product of $f,1_{\mathcal{A}}$ and $\delta$ the sum on the
left hand side splits into $s$ independent sums. We estimate all factors
trivially apart from one of those that contain the balanced function $f$
(which exists by construction of $f_{i}$). Write
$P_{\mathbf{h},\mathbf{j}}=\mathbf{h}+qP_{\mathbf{j}}$ with
$P_{\mathbf{j}}=P_{j_{1}}\times\ldots\times P_{j_{s}}$, then we have for some
$k\leq s$ the inequality

$\displaystyle(\delta P)^{s-1}\Big{|}\sum_{x_{k}\in
h_{k}+qP_{j_{k}}}f(x_{k})\Big{|}\gg\delta^{240s+s+80}N^{s}\gg\delta^{s+80}|P_{j_{k}}|^{s}.$

This implies a density increment of size $\delta\to\delta+\theta\delta^{81}$
for a small $\theta>0$ and a loss of progression length of
$N\to\eta\delta^{240}N$. Now we perform a density increment argument, the
details of which can be found in [3, Section 7].

The number of steps of this iteration is bounded by $\theta^{-1}\delta^{-81}$
and we end up with the condition
$(\eta\delta^{240})^{\theta^{-1}\delta^{-81}}N\ll C_{Q}D^{8}\delta^{-80}$ for
some constant $C_{Q}$, which depends on the smallest non-trivial solution of
$Q$ in $\mathbb{Z}$ (for a discussion of this, see [3, Section 2]). Taking
logarithms and rearranging for $\delta$ gives us $\delta\ll(\log N)^{-1/82}$,
for example.

##  5\. Structure of Quadratic Forms with Off-rank One

Now we turn our attention towards the second topic of this paper, the off-rank
one situation. Given a symmetric matrix with off-rank one, what can we say
about its structure? The following lemma provides a complete answer to this
question. Write $\mathbf{v}\otimes\mathbf{w}:=\mathbf{v}\cdot\mathbf{w}^{T}$
for $\mathbf{v},\mathbf{w}\in\mathbb{R}^{d}$.

######  Lemma 5.1.

Let $Q\in\mathbb{Z}^{s\times s}$ be a symmetric matrix with off-rank $r=1$.
Then either

  * $(i)$

$Q=D+\mathbf{v}\otimes\mathbf{e}_{j}+\mathbf{e}_{j}\otimes\mathbf{v}$, where
$\mathbf{v}\in\mathbb{Z}^{s}$, $\mathbf{e}_{j}$ is the standard basis vector
for some $1\leq j\leq s$ and $D\in\mathbb{Z}^{s\times s}$ is diagonal, or

  * $(ii)$

$Q=m^{-1}(D+\mathbf{v}\otimes\mathbf{v})$ with a diagonal quadric
$D\in\mathbb{Z}^{s\times s}$, some vector $\mathbf{v}\in\mathbb{Z}^{s}$ and
$m\in\mathbb{Z}\backslash\\{0\\}$.

###### Proof.

By permutating variables, if necessary, we can assume that $Q$ has the form

$\displaystyle Q=\begin{pmatrix}g&u&\mathbf{a}^{T}\\\ u&h&\mathbf{b}^{T}\\\
\mathbf{a}&\mathbf{b}&C\\\ \end{pmatrix}$

for some matrix $C\in\mathbb{Z}^{(s-2)\times(s-2)}$ and $u\neq 0$. Now
consider for some $i\neq j$ the off-diagonal submatrix

$\displaystyle\begin{pmatrix}u&a_{j}\\\ b_{i}&c_{ij}\end{pmatrix}.$

For this to have rank at most one, we must have $uc_{ij}=a_{j}b_{i}$. If
$\mathbf{a}=0$ or $\mathbf{b}=0$ this implies that $c_{ij}=0$ for all $i\neq
j$, which corresponds to case (i). Otherwise, we can assume that
$\mathbf{a}\neq 0$ and $\mathbf{b}\neq 0$. Again by the off-rank property we
deduce that the vectors $\mathbf{a}$ and $\mathbf{b}$ are linearly dependent
and we have $\mathbf{b}=\lambda\mathbf{a}$ for some
$\lambda\in\mathbb{Q}\backslash\\{0\\}$. By using $uc_{ij}=a_{j}b_{i}$ we
deduce that $C=u^{-1}\mathbf{b}\otimes\mathbf{a}+E$ for some diagonal matrix
$E$. It remains to ‘lift’ this information to the matrix $Q$.

If we consider the structure of $Q$ modulo diagonal matrices, we have to show
that one can choose diagonal entries $x$ and $y$ in such a way that the matrix

$\displaystyle P=\begin{pmatrix}x&u&\mathbf{a}^{T}\\\
u&y&\lambda\mathbf{a}^{T}\\\
\mathbf{a}&\lambda\mathbf{a}&u^{-1}\lambda\mathbf{a}\otimes\mathbf{a}\\\
\end{pmatrix},$

has the form $m^{-1}\mathbf{v}\otimes\mathbf{v}$ for some
$m\in\mathbb{Z}\backslash\\{0\\}$ and $\mathbf{v}\in\mathbb{Z}^{s}$. There
isn’t any choice but to complete $P$ to

$\displaystyle\begin{pmatrix}u\lambda^{-1}&u&\mathbf{a}^{T}\\\ u&\lambda
u&\lambda\mathbf{a}^{T}\\\
\mathbf{a}&\lambda\mathbf{a}&u^{-1}\lambda\mathbf{a}\otimes\mathbf{a}\\\
\end{pmatrix}$

and check that this matrix has indeed rank one. To be able to write this as a
tensor product of a vector $\mathbf{v}$ with itself, we multiply with the
common factor $m=k^{2}\lambda u$. A suitable $k\in\mathbb{N}$ makes all
entries to integers and allows us to choose $\mathbf{v}^{T}=k\,(u,\lambda
u,\lambda\mathbf{a}^{T})$. ∎

##  6\. The Parameter Method

In this section we are going to analyse case (i) in Lemma 5.1. By changing
notation, we can assume that the ‘non-diagonal’ variable is $x_{s}$ and by
completing the square we can write the quadratic equation $Q(\mathbf{x})=0$ in
the form

$\displaystyle
a_{1}(b_{1}x_{1}-c_{1}x_{s})^{2}+\ldots+a_{s-1}(b_{s-1}x_{s-1}-c_{s-1}x_{s})^{2}+a_{s}x_{s}^{2}=0$

for some $a_{i},b_{i},c_{i}\in\mathbb{Z}$. We can assume that $b_{i}\neq 0$
for all $1\leq i\leq s-1$.

Now we can simplify this equation by incorporating the restriction of
translation invariance. If we replace each $x_{i}$ by $x_{i}+1$ we obtain the
same quadratic form with an additional linear term, which is two times

$\displaystyle
a_{1}(b_{1}x_{1}-c_{1}x_{s})(b_{1}-c_{1})+\ldots+a_{s-1}(b_{s-1}x_{s-1}-c_{s-1}x_{s})(b_{s-1}-c_{s-1})+a_{s}x_{s},$

as well as the constant term

$\displaystyle a_{1}(b_{1}-c_{1})^{2}+\ldots+a_{s-1}(b_{s}-c_{s})^{2}+a_{s}.$

The linear term must be zero for all $\mathbf{x}$, which implies that
$b_{i}=c_{i}$ for all $i$ with $a_{i}\neq 0$. Where $a_{i}=0$ we can simply
assume that $b_{i}=c_{i}$ without changing the quadratic form. Furthermore
this implies that $a_{s}=0$. These conditions are sufficient for the vanishing
of the constant term as well.

If we set $d_{i}=a_{i}b_{i}^{2}$ we end up with the almost diagonal form

(6.1) $\displaystyle
d_{1}(x_{1}-x_{s})^{2}+\ldots+d_{s-1}(x_{s-1}-x_{s})^{2}=0.$

If $d_{i}=0$ for some $i$, we can easily find a non-trivial solution in
$\mathcal{A}$. Therefore, we can assume that all $d_{i}\neq 0$.

We view $x_{s}$ as a free ‘parameter’ ranging over the set $\mathcal{A}$ and
focus on the remaining $s-1$ variables. We define the shifted exponential sum
by

(6.2) $\displaystyle
U_{\mathcal{A}-y}(\alpha)=\sum_{x\in\mathcal{A}}e(\alpha(x-y)^{2}).$

The number of solutions to (6.1) can be written as the Fourier integral

$\displaystyle\sum_{x_{s}\in\mathcal{A}}\int_{0}^{1}\prod_{i=1}^{s-1}U_{\mathcal{A}-x_{s}}(d_{i}\alpha)\,d\alpha.$

We compare this to the situation, where $x_{1},\ldots,x_{s-1}\in[1,N]$. The
number of solutions in this case is given by

$\displaystyle\sum_{x_{s}\in\mathcal{A}}\int_{0}^{1}\prod_{i=1}^{s-1}U_{[1,N]-x_{s}}(d_{i}\alpha)\,d\alpha.$

The inner integral is bounded from below by

(6.3)
$\displaystyle\int_{0}^{1}\prod_{i=1}^{s-1}U_{[1,N/2]}(d_{i}\alpha)\,d\alpha$

since the set $[1,N]-x_{s}$ covers either $[1,N/2]$ or $[-N/2,-1]$ completely
and equation (6.1) has even degree.

As long as $s-1\geq 5$ and not all $d_{i}$ have the same sign, the classical
circle method gives us many integer solutions for a diagonal quadratic
equation. (See the book of Davenport [1], for example.) The case, where all
$d_{i}$ have the same sign is excluded by the assumption that the quadric has
a non-singular real solution. This implies a lower bound of size $N^{s-3}$ for
(6.3) and we obtain

$\displaystyle\sum_{x_{s}\in\mathcal{A}}\int_{0}^{1}\prod_{i=1}^{s-1}U_{[1,N]-x_{s}}(d_{i}\alpha)\,d\alpha\gg\delta
N\cdot N^{s-3}.$

Since there are only $\delta N$ trivial solutions of (6.1) in $\mathcal{A}$ by
assumption, we get

(6.4)
$\displaystyle\sum_{x_{s}\in\mathcal{A}}\int_{0}^{1}\Big{|}\prod_{i=1}^{s-1}U_{\mathcal{A}-x_{s}}(d_{i}\alpha)-\prod_{i=1}^{s-1}\delta
U_{[1,N]-x_{s}}(d_{i}\alpha)\Big{|}\,d\alpha\gg\delta^{s}N^{s-2}-\delta N.$

To deduce a correlation estimate, we replace the indicator function
$1_{\mathcal{A}}$ in (6.2) by $1_{\mathcal{A}}=\delta 1_{[1,N]}+f$ with the
balanced function $f$ from (1.1). An expansion of the first product in (6.4)
creates $2^{s-1}$ terms, where the first one is cancelled by the second
product. The remaining contributions can be bounded from above by a finite sum
of terms of the form

$\displaystyle\sum_{x_{s}\in\mathcal{A}}\int_{0}^{1}\prod_{i=1}^{s-1}|U_{T(x_{s})g_{i}}(d_{i}\alpha)|\,d\alpha$

with $(T_{y}g)(x)=g(x-y)$ and at least one of the $g_{i}$ equal to $f$. By
Hölder’s inequality, we can bound this expression by

$\displaystyle\sup_{x_{s},\alpha}|U_{T_{x_{s}}f}(\alpha)|^{1/2}\sum_{x_{s}\in\mathcal{A}}\int_{0}^{1}|U_{g}(\alpha)|^{s-3/2}\,d\alpha$

for some function $|g|\leq 1$ defined on $[-N,N]$ with $\sum_{|n|\leq
N}|g(n)|\leq 2\delta N$. This implies $|U_{g}(\alpha)|\leq 2\delta N$ and the
upper bound

$\displaystyle(2\delta
N)^{s-6}\sup_{x_{s},\alpha}|U_{T_{x_{s}}f}(\alpha)|^{1/2}\sum_{x_{s}\in\mathcal{A}}\int_{0}^{1}|U_{g}(\alpha)|^{9/2}\,d\alpha.$

We get the bound $O(N^{5/2})$ for the integral from [3, Theorem 6]. This leads
to the estimate

$\displaystyle(\delta
N)^{s-5}\sup_{x_{s},\alpha}|U_{T_{x_{s}}f}(\alpha)|^{1/2}N^{5/2}\gg\delta^{s}N^{s-2}-\delta
N.$

As long as $N\gg_{Q}\delta^{-2}$ we obtain the correlation estimate

$\displaystyle\Big{|}\sum_{x}f(x)e(\alpha(x-x_{s})^{2})\Big{|}\gg\delta^{10}N$

for some $x_{s}\in[1,N]$ and $\alpha\in\mathbb{T}$. By expanding the square,
one can see that the left hand side is just a quadratic exponential sum. The
density increment procedure in [3, Section 7] gives us the bound
$\delta\ll(\log\log N)^{-1/11}$.

##  7\. The Diagonal System

We are given a translation invariant quadratic form
$\mathbf{x}^{T}Q\mathbf{x}=0$ in $s\geq 6$ variables with off-rank one, which
satisfies case $(ii)$ in Lemma 5.1.

We can write the equation $Q(\mathbf{x})=0$ in the form

(7.1) $\displaystyle
d_{1}x_{1}^{2}+\ldots+d_{s}x_{s}^{2}+(v_{1}x_{1}+\ldots+v_{s}x_{s})^{2}=0,$

where the coefficients $d_{i}$ are equal to the diagonal elements in the
matrix $D$ from Lemma 5.1. Consider
$Q(\mathbf{x}+\boldsymbol{1})=Q(\mathbf{x})+2L(\mathbf{x})+Q(\boldsymbol{1})$,
where $L(\mathbf{x})$ is given by

$\displaystyle
d_{1}x_{1}+\ldots+d_{s}x_{s}+(v_{1}x_{1}+\ldots+v_{s}x_{s})(v_{1}+\ldots+v_{s})$

and the constant term $Q(\boldsymbol{1})$ is

$\displaystyle(d_{1}+\ldots+d_{s})+(v_{1}+\ldots+v_{s})^{2}.$

By translation invariance, the linear and constant terms have to disappear for
all $\mathbf{x}\in\mathbb{Z}$. This can only happen if for
$n=v_{1}+\ldots+v_{s}$ we have $d_{i}=-nv_{i}$ and
$d_{1}+\ldots+d_{s}=-n^{2}$. In the special case $n=0$ we end up with the
linear equation

$\displaystyle v_{1}x_{1}+\ldots+v_{s}x_{s}=0,$

which is covered by the method of Roth [4] and gives a better bound than
needed for our theorem here. We assume from now on that $n\neq 0$ and write
equation (7.1) in the form of a system. For a new variable
$h=v_{1}x_{1}+\ldots+v_{s}x_{s}$ with $h\in\mathbb{Z}$ we obtain (after
multiplication of the linear equation by $n$)

(7.2) $\begin{split}d_{1}x_{1}^{2}+\ldots+d_{s}x_{s}^{2}+h^{2}&=0,\\\
d_{1}x_{1}+\ldots+d_{s}x_{s}+nh&=0,\end{split}$

with the condition $d_{1}+\ldots+d_{s}=-n^{2}$. Now we (arbitrarily) restrict
$h$ to $n\mathbb{Z}$, write $h=nx_{0}$ and $d_{0}=n^{2}$. Then system (7.2)
reduces further into

$\begin{split}d_{0}x_{0}^{2}+d_{1}x_{1}^{2}+\ldots+d_{s}x_{s}^{2}=0,\\\
d_{0}x_{0}+d_{1}x_{1}+\ldots+d_{s}x_{s}=0,\end{split}$

with $d_{0}+d_{1}+\ldots+d_{s}=0$.

Now Theorem 2.5 implies the result for the weaker system (7.2). One should
note at this point that conditions $(ii)$ and $(iii)$ for Theorem 2.5 are
implied by the assumptions in Theorem 2.4 or can be assumed to be true without
loss of generality.

##  8\. Remarks and Open Problems

Looking at the cases $r=1$ and $r=5$, where we can deal with equations in
$5+r$ variables, one might have a vague hope to be able to extend this result
to off-ranks $2\leq r\leq 4$. If this were the case, Theorem 1.1 would require
only $s\geq 10$ variables instead of $s\geq 17$. (Recent work of Zhao [6]
could probably reduce the bound even further to $s\geq 9$.)

To prove such a result, we need a structure theorem along the lines of Lemma
5.1 for quadratic forms with off-rank $2\leq r\leq 4$. Let us look at the next
simplest case $r=2$.

When $r=1$ we have two cases. Slightly simplified, they correspond to the
decompositions $M=D+\mathbf{v}\otimes\mathbf{v}$ and
$M=D+\mathbf{e}_{s}\otimes\mathbf{v}+\mathbf{v}\otimes\mathbf{e}_{s}$ for a
diagonal matrix $D$. Combining the two ideas, we end up with three different
structures for $r=2$:

  * (i)

$Q=D+\mathbf{v}\otimes\mathbf{v}+\mathbf{w}\otimes\mathbf{w}$,

  * (ii)

$Q=D+\mathbf{v}\otimes\mathbf{v}+\mathbf{e}_{s}\otimes\mathbf{w}+\mathbf{w}\otimes\mathbf{e}_{s}$,

  * (iii)

$Q=D+\mathbf{e}_{s-1}\otimes\mathbf{v}+\mathbf{v}\otimes\mathbf{e}_{s-1}+\mathbf{e}_{s}\otimes\mathbf{w}+\mathbf{w}\otimes\mathbf{e}_{s}$.

It is easy to check, that in all three cases we have indeed a matrix with off-
rank at most two. Sadly this naïve idea doesn’t work and they don’t cover all
possible cases of matrices with off-rank $r=2$. Consider the matrix

$\displaystyle\begin{pmatrix}*&1&1&0&1\\\ 1&*&0&1&1\\\ 1&0&*&1&0\\\
0&1&1&*&-1\\\ 1&1&0&-1&*\end{pmatrix},$

where stars mark arbitrary entries.

For it to be of type (i), we would need to find diagonal entries such that the
resulting matrix has rank two. By choosing suitable $3\times 3$ matrices with
only one entry missing, we can fill in the diagonal entries easily and check
that this cannot be done consistently.

It is also easy to see that option (iii) is not correct since there are non-
zero entries in more than two rows. For option (ii), we have to show that any
$4\times 4$ submatrix, which results by deleting a row and corresponding
column cannot be completed to have rank one. This follows directly from the
existence of off-diagonal $2\times 2$ matrices of full rank in each of those
cases.

It is an interesting problem, whether this a is pathological counterexample
that can be understood by adding a case (iv) to the above list or whether
symmetric off-rank two matrices don’t have a simple classification.

Even though a complete classification seems a non-trivial task, it is likely,
that the ideas of this paper can be used to improve slightly on the variable
bounds for $r\in\\{2,3,4\\}$ and, therefore, potentially reduce the overall
bound of Theorem 1.1 from $s\geq 17$ to $s\geq 16$, for example.

Concerning the density bounds, it isn’t hard to see that the proof given for
Theorem 2.2 does generalize to any situation in [3], where we use the function
$K(\alpha)$ from (3.1) to bound the $L^{p}$-norm of our exponential sum. This
would take care of almost all quadratic forms with off-rank $1\leq r\leq 4$ as
well. The only cases, where this is not possible correspond to Section 8 in
[3], where we reduce the problem to a linear system in four equations. Further
advances in the linear theory could provide bounds of the form $(\log N)^{-c}$
for all quadratic forms in sufficiently many variables.

## References

  * [1] H. Davenport,  Analytic methods for Diophantine equations and Diophantine inequalities. Second edition.  Cambridge University Press, Cambridge, 2005.
  * [2] E. Keil,  On a diagonal quadric in dense variables, arXiv:1306.4524.
  * [3] E. Keil,  Translation invariant quadratic forms in dense sets, arXiv:1308.6680.
  * [4] K. F. Roth,  On certain sets of integers, J. London Math. Soc. 28, (1953), 104–109.
  * [5] M. Smith,  On solution-free sets for simultaneous quadratic and linear equations, J. Lond. Math. Soc. (2) 79 (2009), no. 2, 273–293.
  * [6] L. Zhao,  The quadratic form in 9 prime variables, arXiv:1402.3697.