# Quantifying the common genetic variability of bacterial traits

T. Tien Mai(1) , Gerry Tonkin-Hill(2),(3) , John A Lees(4),(6)
and Jukka Corander(2),(3),(5)

( (1) Department of Mathematical Sciences, Norwegian University of Science and
Technology, Trondheim, Norway.
(2) Parasites and Microbes, Wellcome Sanger Institute, Cambridgeshire, UK.
(3) Oslo Centre for Biostatistics and Epidemiology, Department of
Biostatistics, University of Oslo, Norway.
(4) MRC Centre for Global Infectious Disease Analysis, Department of
Infectious Disease Epidemiology, Imperial College, UK.
(5) Department of Mathematics and Statistics, University of Helsinki, Finland.
(6)European Molecular Biology Laboratory, European Bioinformatics Institute
EMBL-EBI, Hinxton, UK.
email<EMAIL_ADDRESS>)

###### Abstract

The study of common heritability, or co-heritability, among multiple traits
has been widely established in quantitative and molecular genetics. However,
in bacteria, genome-based estimation of heritability has only been considered
very recently and no methods are currently available for considering co-
heritability. Here we introduce such a method and demonstrate its usefulness
by multi-trait analyses of the three major human pathogens Escherichia coli,
Neisseria gonorrhoeae and Streprococcus pneumoniae. We anticipate that the
increased availability of high-throughput genomic and phenotypic screens of
bacterial populations will spawn ample future opportunities to understand the
common molecular basis of different traits in bacteria.

Keywords: Co-heritability; Heritability; antimicrobial resistance; genetic
relatedness; invasiveness; carriage duration.

##  1 Introduction

Understanding heritability of traits is one of the cornerstones of both
quantitative and molecular genetics, representing a rich history that spans
more than a century of research [Lynch and Walsh, 1998], [Bürger, 2000]. More
recently, the genetic relationships between different traits have received
substantial attention in human molecular genetics, for example leading to the
identification of important relationships such as the one between height and
coronary artery disease, as well as between ulcerative colitis and childhood
obesity [Bulik-Sullivan et al., 2015a]. While methods for genome-wide
association studies (GWAS) originally developed in human genetics have been
successfully adapted to the field of bacterial genomics [Lees et al., 2016,
Collins and Didelot, 2018, Earle et al., 2016], genome-wide estimation of
heritability in bacteria has only more recently received attention [Lees et
al., 2017, Mai et al., 2021]. Thus far, identification of genetic correlations
between traits in bacteria has received little attention. Hence, there is an
opportunity to broaden our understanding of the genetic correlations of
important bacterial traits, including antibiotic resistance, transmission
rate, carriage duration, virulence and gene expression, to name a few.

Compared to studies in human genetics, genetic association studies in bacteria
present unique methodological challenges due in part to the extensive
variation in gene content and sequence diversity even within members of the
same species [Tettelin et al., 2005]. The clonal reproduction mode of bacteria
further implies that linkage disequilibrium (LD) typically stretches much
further than in sexually reproducing species and consequently many genomic
regions may be associated with any given trait [Chen and Shapiro, 2015]. These
issues have been addressed through the use of short reference independent
sequence features as a replacement for the SNP-based analyses common in human
genetics [Lees et al., 2016, Jaillard et al., 2018]. Simulations, linear mixed
models (LMMs) and penalised multiple regression have also been used to account
for the highly clonal nature of bacterial genome datasets [Lees et al., 2016,
Collins and Didelot, 2018, Lees et al., 2020]. Lineage level effects are often
reported when very high LD obscures locus-specific associations [Earle et al.,
2016].

These advances have led to significant discoveries between genetic variation
in bacteria and important traits such as antibiotic resistance [Earle et al.,
2016], carriage duration [Lees et al., 2017] and virulence [Chaguza et al.,
2020, Weinert et al., 2015]. However, investigations of genetic correlations
between such traits have so far mostly been restricted to broader scale
comparisons of lineages, as exemplified by the identification of an
association between antimicrobial resistance and carriage duration in
Streptococcus pneumoniae [Lehtinen et al., 2017].

There exists a large body of literature for investigating genetic correlations
between traits in humans and in animals in general. These include those based
on Mendelian randomisation, which focuses on investigating previous variants
found to be significantly associated with a phenotype, to investigate causal
relationships between risk factors and disease. Other approaches include
restricted maximum likelihood (REML) and polygenic scores which typically make
use of linear mixed models (LMMs) to account for all observed SNPs. A
comprehensive discussion of these approaches can be found in Van Rheenen et
al., [van Rheenen et al., 2019]. For many complex human traits such as height,
heritability is spread over many sites with each contributing a small effect
[Yang et al., 2010, Speed et al., 2017, Speed and Balding, 2019]. Conversely,
in bacteria, many phenotypes of interest are the result of a small number of
variants with larger effect sizes, with a notable example being the resistance
conferring mutations in the penicillin binding proteins of S. pneumoniae
[Chewapreecha et al., 2014b]. For this reason, in bacterial genomics,
penalised regression techniques present an attractive alternative to the more
common LMM approach as it has been shown to better account for the clonal
population structure of bacteria [Saber and Shapiro, 2020, Lees et al., 2020].

Here, we develop a high-dimensional penalised linear regression model to
investigate genetic correlations between bacterial phenotypes. We demonstrate
the accuracy and advantages of the approach by investigating the association
between resistance to different antibiotic classes, carriage duration and
virulence in a selection of Streptococcus pneumoniae, Escherichia coli and
Neisseria gonorrhoeae data sets.

##  2 Results

###  2.1 Overview

To investigate the genetic covariance and coheritability between two
phenotypic traits of interest, we built off of previous work which used
elastic net regression to generate accurate phenotypic prediction models and
more accurately estimated heritability in related populations [Lees et al.,
2020]. Here, we model each phenotype as a linear combination of genetic
variants with independent error terms to account for environmental and
unmeasured genetic effects. The genetic covariance between the two traits is
then obtained by subtracting the covariance of the inferred error terms from
the covariance between the two traits of interest (Methods). This can then be
standardised to obtain the genetic correlation with the resulting scale being
between -1 indicating perfect negative correlation and 1 indicating perfect
positive correlation between the estimated genetic effects. Similar to
previous studies in human genetics, in addition to genetic correlation, we
also consider the covariance of the causal effects [Shi et al., 2017, Bulik-
Sullivan et al., 2015b]. This allows for the separation of causal correlations
from those that may be driven by linkage disequilibrium (LD) and is
particularly important in analysing highly clonal datasets such as those
commonly found in bacterial genomics. Using this approach, we examined the
genetic correlations between resistance to different antibiotics across a
number of bacterial species as well as virulence and carriage duration in
Streptococcus pneumoniae.

###  2.2 Antibiotic resistance, virulence and carriage duration are
genetically correlated in Streptococcus pneumoniae

To examine the ability of our approach to identify meaningful genetic
correlations between bacterial phenotypes, we initially considered two large
S. pneumoniae datasets. The Maela dataset consists of 3069 whole genomes taken
from an infant cohort study in a refugee camp on the Thailand-Myanmar border
[Chewapreecha et al., 2014a]. This dataset has been used in genome wide
association studies to identify genetic loci associated with antibiotic
resistance pneumococcal carriage duration [Lees et al., 2017, Chewapreecha et
al., 2014b]. Using our approach, we investigated the genetic correlation
between the resistance to five antibiotics and carriage duration. Figure 1
indicates that pneumococcal carriage duration is positively correlated with
resistance particularly in the case of Erythromycin. We found a similar
correlation in the cases of Penicillin, Tetracycline and co-trimoxazole
although there was a higher degree of uncertainty in these estimates. This is
consistent with a recent study that investigated the association between
resistance and carriage duration in the Maela dataset [Lehtinen et al., 2017].
Interestingly, we found that there appears to be no covariance between the
causal effects of carriage duration and resistance to any of the antibiotics
considered. While laboratory studies have shown that resistance to antibiotics
leads to reduced growth rates in S. pneumoniae [Rozen et al., 2007, Trzciński
et al., 2006] these results suggest that it is not the fitness costs of
resistance that drive the association with carriage duration. Rather the
fitness benefits of resistance are likely to be greater for those lineages
with a longer duration of carriage due to the increased probability that they
encounter antimicrobial treatment.

In addition to the correlation with carriage duration, we found that there
were strong genetic correlations between resistance to each of the antibiotics
considered (Figure 1). In the cases of tetracycline resistance, erythromycin
resistance, and chloramphenicol resistance, we also observed correlations
between the causal effects. To investigate whether these correlations were
unique to this dataset, we ran the same analysis on a large subset of the
Global Pneumococcal Sequencing project consisting of genomes sampled in South
Africa for which accurate resistance and virulence phenotypic information was
known [Lo et al., 2019, Gladstone et al., 2019]. Figure 2 shows that a very
similar profile of genetic correlations between resistance to the same
antibiotics was present in this dataset, which indicates that our results are
consistent across multiple distinct locations. Although resistance to these
antibiotics is thought to be mediated by different mechanisms, the
corresponding resistance genes can be found in a single cassette ICESp23FST81.
This cassette can be inherited both vertically and horizontally and varies in
which genes it harbours: tetM for tetracycline; ermB or mel and mef for
erythromycin; cat for chloramphenicol [Croucher et al., 2011]. Thus, the
causal correlation identified is likely to be driven by the presence and
absence of this cassette which is not in strong LD with the rest of the
genome.

Figure 1: Maela data. Heritability in diagonals, genetic relatedness in the
upper diagonals, genetic correlation in lower diagonals.

The South African dataset also included accurate information on which genomes
were isolated from invasive disease cases allowing us to consider the genetic
correlation between resistance and virulence. We found the virulence was
negatively genetically correlated with resistance although there was
considerable uncertainty in these estimates. Similar to resistance and
carriage duration, we did not observe a strong covariance between the causal
effects of virulence and resistance. The anti-correlation between virulence
and resistance has been observed previously in isolates of S. pneumoniae from
China where resistance was observed at much lower frequencies in isolates
taken from children with meningitis [Wang et al., 2019]. A similar finding has
also been made in Klebsiella pneumoniae where common virulence genes are
rarely found together with antimicrobial resistance genes [Holt et al., 2015].
A potential explanation for this is that, in contrast to carriage duration,
resistance is less beneficial to those pneumococcal lineages that commonly
cause invasive disease. In S. pneumoniae, these lineages are often rarely
found in carriage studies indicating they may encounter antimicrobial
treatment less frequently than the common multidrug resistant lineages
[Chaguza et al., 2016, Gladstone et al., 2019].

Figure 2: Streptococcus pneumoniae data. Heritability in diagonals, genetic
relatedness in the upper diagonals, genetic correlation in lower diagonals.

###  2.3 Genetic correlation between resistance to different antibiotics is
observed in multiple bacterial species

We also observed genetic correlations between resistance phenotypes in both
Escherichia coli and Neisseria gonorrhoeae. Using a dataset of 1509
Escherichia coli isolates taken from an 11 year systematic hospital survey of
bacteremia associated isolates in England, we explored whether there were
identifiable correlations between the resistance profiles of each isolate
[Kallonen et al., 2017]. Figure 3 demonstrates that with the exception of
Amoxicillin there was a positive genetic correlation between the resistance
phenotypes. However, we observed considerable uncertainty in these estimates
and thus can not rule out correlations between the other resistance profiles.

We also investigated the genetic correlations between antibiotic resistance
phenotypes in 1595 N. gonorrhoeae isolates collected from the USA, Canada and
England [Unemo et al., 2016, Grad et al., 2016, Demczuk et al., 2015, De Silva
et al., 2016, Schubert et al., 2018]. This identified strong positive genetic
correlations between the resistance phenotypes to Tetracycline, Ciproflaxcin
and Penicillin (Figure 4). Similar to S. pneumoniae, N. gonorrhoeae is thought
to have different resistance mechanisms to each of these antibiotics. However,
unlike in S. pneumoniae we did not find a causal association between the
effects of each antibiotic. Rather it is likely that the genetic correlation
is driven by the two major modern gonococcal lineages with a multidrug
resistance lineage being common in high-risk sexual networks and a
multisusceptible lineage often associated with lower risk heterosexual
networks [Sánchez-Busó et al., 2019]. Strikingly, we also observed a negative
genetic correlation between azithromycin and tetracycline or ciprofloxacin
albeit with a greater level of uncertainty. This is consistent with the very
rare instances of resistance to the common dual therapy (injectable
ceftriaxone plus oral azithromycin) suggesting that the fitness costs of
gaining resistance to azithromycin in addition to certain other antibiotics is
sufficiently high to prevent such isolates from proliferating [Sánchez-Busó et
al., 2019, Grad et al., 2016].

Figure 3: Ecoli data. Heritability in diagonals, genetic covariance in the
upper diagonals, genetic correlation in lower diagonals. Figure 4: Neisseria
gonorrhoeae data. Heritability in diagonals, genetic covariance in the upper
diagonals, genetic correlation in lower diagonals.

##  3 Discussion

The increasingly common use of genome wide association studies in bacterial
genomics has led to vast improvements in our understanding of the links
between bacterial genomes and important phenotypes such as antibiotic
resistance and virulence. Here, we have developed a method to estimate the
coheritability between traits in bacterial populations. Building on our
earlier work using elastic net regression to fit whole genome models to
estimate regression slopes [Lees et al., 2020], and our approaches to more
accurately estimate heritability in related populations [Mai et al., 2021], we
are able to use the covariance between fitted predictor values to estimate
shared genetic contributions between pairs of traits. By looking at causal
correlations, we have demonstrated that it is possible to distinguish when two
phenotypes are likely to be the result of similar genetic mutations or when
population structure and epidemiology could be driving the genetic
correlation. This is highlighted by the observed associations between
resistance, virulence and carriage duration in S. pneumoniae. Here,
statistical genetics approaches are particularly powerful as these phenotypes
can be difficult to measure in lab models, whereas they can be more readily
quantified in the natural population. We show that consistent with their known
causal genetic mechanisms, typically point mutations or genes which are not
co-localised, these phenotypes are unlikely to be driven by the same genetic
features.

In S. pneumoniae we found evidence of co-heritability between tetracycline
resistance, erythromycin resistance, and chloramphenicol resistance – though
the strength and significance of this finding varied between populations.
Although the causal elements themselves are not directly coheritable, the
molecule they can be found on is inherited both vertically and horizontally.
From a more methodological perspective, the SNPs on these elements would
typically all be associated with each phenotype, and shared when multiple
resistances are present [Lees et al., 2016]. The replication of results across
distinct pneumococcal datasets in Thailand and South Africa suggests that the
method is robust to the impacts of population structure and sampling biases.
We also demonstrated weaker genetic correlations between resistance to
antibiotics in N. gonorrhoeae and E. coli. Unlike in S. pneumoniae, these were
found to not have causal correlations and are likely to be the result of the
significant resistance associated with population structure observed in the
species.

In summary, our method is capable of quantifying the shared genetic
variability related to variance in pairs of phenotypes, in closely related
bacterial populations. We did not explicitly model LD or relatedness, nor look
at accessory genome variation or more complex genetic variation. Additionally,
our findings at this point are largely observational. Nevertheless, we believe
this will be a useful first step into further research on shared genomic bases
between traits in pathogenic bacteria. To enable other researchers to easily
conduct similar analyses, we have also provided an accompanying R package on
GitHub (https://github.com/tienmt/coher).

## Acknowledgments

TTM is supported by the Norwegian Research Council grant number 309960 through
the Centre for Geophysical Forecasting at NTNU. JAL acknowledges funding from
the MRC Centre for Global Infectious Disease Analysis (reference
MR/R015600/1), jointly funded by the UK Medical Research Council (MRC) and the
UK Foreign, Commonwealth & Development Office (FCDO), under the MRC/FCDO
Concordat agreement and is also part of the EDCTP2 programme supported by the
European Union.

## References

  * [Belloni et al., 2013] Belloni, A., Chernozhukov, V., et al. (2013).  Least squares after model selection in high-dimensional sparse models.  Bernoulli, 19(2):521–547.
  * [Bulik-Sullivan et al., 2015a] Bulik-Sullivan, B., Finucane, H. K., Anttila, V., Gusev, A., Day, F. R., Loh, P.-R., ReproGen Consortium, Psychiatric Genomics Consortium, Genetic Consortium for Anorexia Nervosa of the Wellcome Trust Case Control Consortium 3, Duncan, L., Perry, J. R. B., Patterson, N., Robinson, E. B., Daly, M. J., Price, A. L., and Neale, B. M. (2015a).  An atlas of genetic correlations across human diseases and traits.  Nat. Genet., 47(11):1236–1241.
  * [Bulik-Sullivan et al., 2015b] Bulik-Sullivan, B. K., Loh, P.-R., Finucane, H. K., Ripke, S., Yang, J., Patterson, N., Daly, M. J., Price, A. L., Neale, B. M., of the Psychiatric Genomics Consortium, S. W. G., et al. (2015b).  Ld score regression distinguishes confounding from polygenicity in genome-wide association studies.  Nature genetics, 47(3):291.
  * [Bürger, 2000] Bürger, R. (2000).  The mathematical theory of selection, recombination, and mutation.  John Wiley & Sons.
  * [Chaguza et al., 2016] Chaguza, C., Cornick, J. E., Harris, S. R., Andam, C. P., Bricio-Moreno, L., Yang, M., Yalcin, F., Ousmane, S., Govindpersad, S., Senghore, M., Ebruke, C., Du Plessis, M., Kiran, A. M., Pluschke, G., Sigauque, B., McGee, L., Klugman, K. P., Turner, P., Corander, J., Parkhill, J., Collard, J.-M., Antonio, M., von Gottberg, A., Heyderman, R. S., French, N., Kadioglu, A., Hanage, W. P., Everett, D. B., Bentley, S. D., and PAGe Consortium (2016).  Understanding pneumococcal serotype 1 biology through population genomic analysis.  BMC Infect. Dis., 16(1):649.
  * [Chaguza et al., 2020] Chaguza, C., Yang, M., Cornick, J. E., du Plessis, M., Gladstone, R. A., Kwambana-Adams, B. A., Lo, S. W., Ebruke, C., Tonkin-Hill, G., Peno, C., Senghore, M., Obaro, S. K., Ousmane, S., Pluschke, G., Collard, J.-M., Sigaùque, B., French, N., Klugman, K. P., Heyderman, R. S., McGee, L., Antonio, M., Breiman, R. F., von Gottberg, A., Everett, D. B., Kadioglu, A., and Bentley, S. D. (2020).  Bacterial genome-wide association study of hyper-virulent pneumococcal serotype 1 identifies genetic variation associated with neurotropism.  Commun Biol, 3(1):559.
  * [Chen and Shapiro, 2015] Chen, P. E. and Shapiro, B. J. (2015).  The advent of genome-wide association studies for bacteria.  Curr. Opin. Microbiol., 25:17–24.
  * [Chewapreecha et al., 2014a] Chewapreecha, C., Harris, S. R., Croucher, N. J., Turner, C., Marttinen, P., Cheng, L., Pessia, A., Aanensen, D. M., Mather, A. E., Page, A. J., Salter, S. J., Harris, D., Nosten, F., Goldblatt, D., Corander, J., Parkhill, J., Turner, P., and Bentley, S. D. (2014a).  Dense genomic sampling identifies highways of pneumococcal recombination.  Nat. Genet., 46(3):305–309.
  * [Chewapreecha et al., 2014b] Chewapreecha, C., Marttinen, P., Croucher, N. J., Salter, S. J., Harris, S. R., Mather, A. E., Hanage, W. P., Goldblatt, D., Nosten, F. H., Turner, C., Turner, P., Bentley, S. D., and Parkhill, J. (2014b).  Comprehensive identification of single nucleotide polymorphisms associated with beta-lactam resistance within pneumococcal mosaic genes.  PLoS Genet., 10(8):e1004547.
  * [Collins and Didelot, 2018] Collins, C. and Didelot, X. (2018).  A phylogenetic method to perform genome-wide association studies in microbes that accounts for population structure and recombination.  PLoS Comput. Biol., 14(2):e1005958.
  * [Croucher et al., 2011] Croucher, N. J., Harris, S. R., Fraser, C., Quail, M. A., Burton, J., van der Linden, M., McGee, L., von Gottberg, A., Song, J. H., Ko, K. S., Pichon, B., Baker, S., Parry, C. M., Lambertsen, L. M., Shahinas, D., Pillai, D. R., Mitchell, T. J., Dougan, G., Tomasz, A., Klugman, K. P., Parkhill, J., Hanage, W. P., and Bentley, S. D. (2011).  Rapid pneumococcal evolution in response to clinical interventions.  Science, 331(6016):430–434.
  * [De Silva et al., 2016] De Silva, D., Peters, J., Cole, K., Cole, M. J., Cresswell, F., Dean, G., Dave, J., Thomas, D. R., Foster, K., Waldram, A., Wilson, D. J., Didelot, X., Grad, Y. H., Crook, D. W., Peto, T. E., Walker, A. S., Paul, J., and Eyre, D. W. (2016).  Whole-genome sequencing to determine transmission of neisseria gonorrhoeae: an observational study.  Lancet Infect. Dis., 16(11):1295–1303.
  * [Demczuk et al., 2015] Demczuk, W., Lynch, T., Martin, I., Van Domselaar, G., Graham, M., Bharat, A., Allen, V., Hoang, L., Lefebvre, B., Tyrrell, G., Horsman, G., Haldane, D., Garceau, R., Wylie, J., Wong, T., and Mulvey, M. R. (2015).  Whole-genome phylogenomic heterogeneity of neisseria gonorrhoeae isolates with decreased cephalosporin susceptibility collected in canada between 1989 and 2013.  J. Clin. Microbiol., 53(1):191–200.
  * [Earle et al., 2016] Earle, S. G., Wu, C.-H., Charlesworth, J., Stoesser, N., Gordon, N. C., Walker, T. M., Spencer, C. C. A., Iqbal, Z., Clifton, D. A., Hopkins, K. L., Woodford, N., Smith, E. G., Ismail, N., Llewelyn, M. J., Peto, T. E., Crook, D. W., McVean, G., Walker, A. S., and Wilson, D. J. (2016).  Identifying lineage effects when controlling for population structure improves power in bacterial association studies.  Nat Microbiol, 1:16041.
  * [Gladstone et al., 2019] Gladstone, R. A., Lo, S. W., Lees, J. A., Croucher, N. J., van Tonder, A. J., Corander, J., Page, A. J., Marttinen, P., Bentley, L. J., Ochoa, T. J., Ho, P. L., du Plessis, M., Cornick, J. E., Kwambana-Adams, B., Benisty, R., Nzenze, S. A., Madhi, S. A., Hawkins, P. A., Everett, D. B., Antonio, M., Dagan, R., Klugman, K. P., von Gottberg, A., McGee, L., Breiman, R. F., and Bentley, S. D. (2019).  International genomic definition of pneumococcal lineages, to contextualise disease, antibiotic resistance and vaccine impact.  EBioMedicine, 43:338–346.
  * [Grad et al., 2016] Grad, Y. H., Harris, S. R., Kirkcaldy, R. D., Green, A. G., Marks, D. S., Bentley, S. D., Trees, D., and Lipsitch, M. (2016).  Genomic epidemiology of gonococcal resistance to Extended-Spectrum cephalosporins, macrolides, and fluoroquinolones in the united states, 2000–2013.  J. Infect. Dis., 214(10):1579–1587.
  * [Holt et al., 2015] Holt, K. E., Wertheim, H., Zadoks, R. N., Baker, S., Whitehouse, C. A., Dance, D., Jenney, A., Connor, T. R., Hsu, L. Y., Severin, J., Brisse, S., Cao, H., Wilksch, J., Gorrie, C., Schultz, M. B., Edwards, D. J., Van Nguyen, K., Nguyen, T. V., Dao, T. T., Mensink, M., Minh, V. L., Nhu, N. T. K., Schultsz, C., Kuntaman, K., Newton, P. N., Moore, C. E., Strugnell, R. A., and Thomson, N. R. (2015).  Genomic analysis of diversity, population structure, virulence, and antimicrobial resistance in klebsiella pneumoniae, an urgent threat to public health.  Proc. Natl. Acad. Sci. U. S. A., 112(27):E3574–81.
  * [Jaillard et al., 2018] Jaillard, M., Lima, L., Tournoud, M., Mahé, P., van Belkum, A., Lacroix, V., and Jacob, L. (2018).  A fast and agnostic method for bacterial genome-wide association studies: Bridging the gap between k-mers and genetic events.  PLoS Genet., 14(11):e1007758.
  * [Kallonen et al., 2017] Kallonen, T., Brodrick, H. J., Harris, S. R., Corander, J., Brown, N. M., Martin, V., Peacock, S. J., and Parkhill, J. (2017).  Systematic longitudinal survey of invasive escherichia coli in england demonstrates a stable population structure only transiently disturbed by the emergence of ST131.  Genome Res.
  * [Lees et al., 2017] Lees, J. A., Croucher, N. J., Goldblatt, D., Nosten, F., Parkhill, J., Turner, C., Turner, P., and Bentley, S. D. (2017).  Genome-wide identification of lineage and locus specific variation associated with pneumococcal carriage duration.  Elife, 6.
  * [Lees et al., 2018] Lees, J. A., Galardini, M., Bentley, S. D., Weiser, J. N., and Corander, J. (2018).  pyseer: a comprehensive tool for microbial pangenome-wide association studies.  Bioinformatics, 34(24):4310–4312.
  * [Lees et al., 2020] Lees, J. A., Mai, T. T., Galardini, M., Wheeler, N. E., Horsfield, S. T., Parkhill, J., and Corander, J. (2020).  Improved prediction of bacterial Genotype-Phenotype associations using interpretable Pangenome-Spanning regressions.  MBio, 11(4).
  * [Lees et al., 2016] Lees, J. A., Vehkala, M., Välimäki, N., Harris, S. R., Chewapreecha, C., Croucher, N. J., Marttinen, P., Davies, M. R., Steer, A. C., Tong, S. Y. C., Honkela, A., Parkhill, J., Bentley, S. D., and Corander, J. (2016).  Sequence element enrichment analysis to determine the genetic basis of bacterial phenotypes.  Nat. Commun., 7:12797.
  * [Lehtinen et al., 2017] Lehtinen, S., Blanquart, F., Croucher, N. J., Turner, P., Lipsitch, M., and Fraser, C. (2017).  Evolution of antibiotic resistance is linked to any genetic mechanism affecting bacterial duration of carriage.  Proc. Natl. Acad. Sci. U. S. A., 114(5):1075–1080.
  * [Lo et al., 2019] Lo, S. W., Gladstone, R. A., van Tonder, A. J., Lees, J. A., du Plessis, M., Benisty, R., Givon-Lavi, N., Hawkins, P. A., Cornick, J. E., Kwambana-Adams, B., Law, P. Y., Ho, P. L., Antonio, M., Everett, D. B., Dagan, R., von Gottberg, A., Klugman, K. P., McGee, L., Breiman, R. F., Bentley, S. D., and Global Pneumococcal Sequencing Consortium (2019).  Pneumococcal lineages associated with serotype replacement and antibiotic resistance in childhood invasive pneumococcal disease in the post-PCV13 era: an international whole-genome sequencing study.  Lancet Infect. Dis., 19(7):759–769.
  * [Lynch and Walsh, 1998] Lynch, M. and Walsh, B. (1998).  Genetics and analysis of quantitative traits, volume 1.  Sinauer Sunderland, MA.
  * [Mai et al., 2021] Mai, T. T., Turner, P., and Corander, J. (2021).  Boosting heritability: estimating the genetic component of phenotypic variation with multiple sample splitting.  BMC Bioinformatics, 22(1):164.
  * [Rozen et al., 2007] Rozen, D. E., McGee, L., Levin, B. R., and Klugman, K. P. (2007).  Fitness costs of fluoroquinolone resistance in streptococcus pneumoniae.
  * [Saber and Shapiro, 2020] Saber, M. M. and Shapiro, B. J. (2020).  Benchmarking bacterial genome-wide association study methods using simulated genomes and phenotypes.  Microb Genom, 6(3).
  * [Sánchez-Busó et al., 2019] Sánchez-Busó, L., Golparian, D., Corander, J., Grad, Y. H., Ohnishi, M., Flemming, R., Parkhill, J., Bentley, S. D., Unemo, M., and Harris, S. R. (2019).  The impact of antimicrobials on gonococcal evolution.  Nat Microbiol, 4(11):1941–1950.
  * [Schubert et al., 2018] Schubert, B., Maddamsetti, R., Nyman, J., Farhat, M. R., and Marks, D. S. (2018).  Genome-wide discovery of epistatic loci affecting antibiotic resistance in neisseria gonorrhoeae using evolutionary couplings.  Nature Microbiology, 4(2):328–338.
  * [Shi et al., 2017] Shi, H., Mancuso, N., Spendlove, S., and Pasaniuc, B. (2017).  Local genetic correlation gives insights into the shared genetic architecture of complex traits.  Am. J. Hum. Genet., 101(5):737–751.
  * [Speed and Balding, 2019] Speed, D. and Balding, D. J. (2019).  Sumher better estimates the snp heritability of complex traits from summary statistics.  Nature genetics, 51(2):277.
  * [Speed et al., 2017] Speed, D., Cai, N., Johnson, M. R., Nejentsev, S., Balding, D. J., Consortium, U., et al. (2017).  Reevaluation of snp heritability in complex human traits.  Nature genetics, 49(7):986.
  * [Tettelin et al., 2005] Tettelin, H., Masignani, V., Cieslewicz, M. J., Donati, C., Medini, D., Ward, N. L., Angiuoli, S. V., Crabtree, J., Jones, A. L., Durkin, A. S., Deboy, R. T., Davidsen, T. M., Mora, M., Scarselli, M., Margarit y Ros, I., Peterson, J. D., Hauser, C. R., Sundaram, J. P., Nelson, W. C., Madupu, R., Brinkac, L. M., Dodson, R. J., Rosovitz, M. J., Sullivan, S. A., Daugherty, S. C., Haft, D. H., Selengut, J., Gwinn, M. L., Zhou, L., Zafar, N., Khouri, H., Radune, D., Dimitrov, G., Watkins, K., O’Connor, K. J. B., Smith, S., Utterback, T. R., White, O., Rubens, C. E., Grandi, G., Madoff, L. C., Kasper, D. L., Telford, J. L., Wessels, M. R., Rappuoli, R., and Fraser, C. M. (2005).  Genome analysis of multiple pathogenic isolates of streptococcus agalactiae: implications for the microbial “pan-genome”.  Proc. Natl. Acad. Sci. U. S. A., 102(39):13950–13955.
  * [Trzciński et al., 2006] Trzciński, K., Thompson, C. M., Gilbey, A. M., Dowson, C. G., and Lipsitch, M. (2006).  Incremental increase in fitness cost with increased $\beta$-Lactam resistance in pneumococci evaluated by competition in an infant rat nasal colonization model.  J. Infect. Dis., 193(9):1296–1303.
  * [Unemo et al., 2016] Unemo, M., Golparian, D., Sánchez-Busó, L., and others (2016).  The novel 2016 WHO neisseria gonorrhoeae reference strains for global quality assurance of laboratory investigations: phenotypic, genetic and reference genome ….  Journal of.
  * [van Rheenen et al., 2019] van Rheenen, W., Peyrot, W. J., Schork, A. J., Lee, S. H., and Wray, N. R. (2019).  Genetic correlations of polygenic disease traits: from theory to practice.  Nat. Rev. Genet., 20(10):567–581.
  * [Wang et al., 2019] Wang, C.-Y., Chen, Y.-H., Fang, C., Zhou, M.-M., Xu, H.-M., Jing, C.-M., Deng, H.-L., Cai, H.-J., Jia, K., Han, S.-Z., Yu, H., Wang, A.-M., Yin, D.-D., Wang, C.-Q., Wang, W., Huang, W.-C., Deng, J.-K., Zhao, R.-Z., Chen, Y.-P., Yang, J.-H., Wang, C., Che, Y.-R., Nie, X.-Z., Wang, S.-F., Hao, J.-H., and Zhang, C.-H. (2019).  Antibiotic resistance profiles and multidrug resistance patterns of streptococcus pneumoniae in pediatrics: A multicenter retrospective study in mainland china.  Medicine, 98(24):e15942.
  * [Weinert et al., 2015] Weinert, L. A., Chaudhuri, R. R., Wang, J., Peters, S. E., Corander, J., Jombart, T., Baig, A., Howell, K. J., Vehkala, M., Välimäki, N., Harris, D., Chieu, T. T. B., Van Vinh Chau, N., Campbell, J., Schultsz, C., Parkhill, J., Bentley, S. D., Langford, P. R., Rycroft, A. N., Wren, B. W., Farrar, J., Baker, S., Hoa, N. T., Holden, M. T. G., Tucker, A. W., Maskell, D. J., and BRaDP1T Consortium (2015).  Genomic signatures of human and animal disease in the zoonotic pathogen streptococcus suis.  Nat. Commun., 6:6740.
  * [Yang et al., 2010] Yang, J., Benyamin, B., McEvoy, B. P., Gordon, S., Henders, A. K., Nyholt, D. R., Madden, P. A., Heath, A. C., Martin, N. G., Montgomery, G. W., Goddard, M. E., and Visscher, P. M. (2010).  Common SNPs explain a large proportion of the heritability for human height.  Nat. Genet., 42(7):565–569.
  * [Zou and Hastie, 2005] Zou, H. and Hastie, T. (2005).  Regularization and variable selection via the elastic net.  Journal of the royal statistical society: series B (statistical methodology), 67(2):301–320.

##  Appendix A Methods

###  A.1 Model and definitions

Given two traits $y$ and $z$, each is modelled as a linear combination of $p$
genetic variants $X_{.j}$ and an error term (environmental and unmeasured
genetic effects)

$\displaystyle y_{n}=X_{n\times p}\beta_{p}+\varepsilon_{n},$ (1)
$\displaystyle z_{n}=X_{n\times p}\alpha_{p}+\gamma_{n}.$

We assume that $X_{i.}$ are i.i.d random vector with 0-mean and covariance
matrix $\Sigma$ and that the random noises $\varepsilon_{i},\gamma_{i}$ are
independent of $X$ with 0-mean and the variances
$\sigma_{\varepsilon}^{2},\sigma_{\gamma}^{2}$. We assume that $y$ and $z$
have 0-mean.

######  Remark.

In this work, we do not assume the effects $\beta,\alpha$ follow any
distribution. This is different with most works in the literature that
conducted in linear mixed model. In linear mixed model, the elements of
$\beta$ (and $\alpha$) are considered as i.i.d random variables following a
Gaussian distribution i.e
$\beta_{j}\overset{i.i.d}{\sim}\mathcal{N}(0,\sigma^{2}_{\beta})$, while the
genetic covariates $X$ are assumed fixed.

#### Genetic covariance and correlation

From the problem formulation (1), the covariance between the phenotypes $y$
and $z$ can be explained as a summation of the genetic covariance and
environmental covariance as follow

$\displaystyle{\rm Cov}(y,z)$
$\displaystyle=\mathbb{E}(yz)-\mathbb{E}(y)\mathbb{E}(z)$
$\displaystyle=\mathbb{E}\left\langle
X\beta+\varepsilon,X\alpha+\gamma\right\rangle$
$\displaystyle=\alpha^{\top}\mathbb{E}\left(X^{\top}X\right)\beta+\mathbb{E}(\varepsilon\gamma)$
$\displaystyle=\alpha^{\top}\Sigma\beta+{\rm Cov}(\varepsilon,\gamma).$

As a consequence, the genetic covariance between two traits $y$ and $z$ is
defined by

$\displaystyle\boxed{{\rm g.Cov}(y,z)=\alpha^{\top}\Sigma\beta},$ (2)

and the genetic correlation of two traits $y$ and $z$ is defined by

$\displaystyle\boxed{{\rm
g.Cor}(y,z)=\frac{\alpha^{\top}\Sigma\beta}{\sqrt{\alpha^{\top}\Sigma\alpha}\sqrt{\beta^{\top}\Sigma\beta}}},$
(3)

which is a standardization of the genetic covariance.

The quantities $\beta^{\top}\Sigma\beta$ and $\alpha^{\top}\Sigma\alpha$ are
the genetic variance of each trait and they are used to define the
heritability of the trait.

######  Remark.

Our formula in the definition of the genetic correlation (3) is different to
the one defined in [Shi et al., 2017] where they define the denominator as the
co-heritability . Our formula guarantees that the genetic correlation varies
between -1 and 1 whereas [Shi et al., 2017] have to threshold it.

#### Heritability

Heritability is an important quantitative concept that measures the variation
observed in a trait contributed by the variation of genotype. From (1), we
have

$\displaystyle{\rm
Var}(y_{i})=\beta^{\top}\Sigma\beta+\sigma_{\varepsilon}^{2};\quad{\rm
Var}(z_{i})=\alpha^{\top}\Sigma\alpha+\sigma_{\gamma}^{2}.$

The (narrow-sense) heritabilities of $y$ and $z$ are defined as

$\displaystyle\boxed{h^{2}(y)=\frac{\beta^{\top}\Sigma\beta}{\beta^{\top}\Sigma\beta+\sigma_{\varepsilon}^{2}};\quad
h^{2}(z)=\frac{\alpha^{\top}\Sigma\alpha}{\alpha^{\top}\Sigma\alpha+\sigma_{\gamma}^{2}}.}$
(4)

#### Genetic relatedness

The genetic relatedness as in [Bulik-Sullivan et al., 2015a], or the
covariance of the causal effect as in [Shi et al., 2017], between $y$ and $z$
is defined as the inner product of the regression coefficients

$\displaystyle\boxed{{\rm g.R}(y,z)=\left\langle\beta,\alpha\right\rangle}.$
(5)

We further consider the normalization of the genetic relatedness as it helps
to compare across the phenotypes as well as across the region, it is defined
as

$\boxed{{\rm
g.R}(y,z)=\frac{\left\langle\beta,\alpha\right\rangle}{\|\alpha\|\|\beta\|}}.$
(6)

######  Remark.

Note that in the special case where there is no LD (i.e. the covariance matrix
$\Sigma$ is the identity matrix), the genetic relatedness is exactly the
genetic covariance: ${\rm
g.Cov}(y,z)=\beta^{\top}\Sigma\alpha=\beta^{\top}\mathbf{I}\alpha=\left\langle\beta,\alpha\right\rangle={\rm
g.R}(y,z)$. To make it more clear, consider a toy example with only 2 genetic
variants, let $\beta=(1,0)^{\top}$ and $\alpha=(0,1)^{\top}$ be the causal
effects of the 2 phenotypes $y$ and $z$ and

$\Sigma=\begin{bmatrix}1&0.8\\\ 0.8&1\end{bmatrix}.$

We obtain ${\rm g.R}(y,z)=\left\langle\beta,\alpha\right\rangle=0$ while ${\rm
g.Cov}(y,z)=\beta^{\top}\Sigma\alpha=0.8$. Therefore, when the causal variants
are completely distinct for 2 traits, the genetic relatedness is always zero,
while the genetic covariance can be different from zero due to the LD. This
has also been discussed in [Shi et al., 2017].

###  A.2 Estimation methods

A naive approach to calculate genetic covariance and correlation is to first
estimate the effect sizes $\beta,\alpha$ and the covariance matrix $\Sigma$
then plug-in these values into their formula.

#### Estimating the effect sizes

To estimate the regression coefficients $\beta$ (and $\alpha$), we use the
Elastic net [Zou and Hastie, 2005]:

$\displaystyle\hat{\beta}_{enet}:=\arg\min_{\beta_{0},\beta}\frac{1}{n}\sum_{i=1}^{n}\ell(y_{i},\beta_{0}+\beta^{T}x_{i})+\lambda\left[0.5(1-\alpha)||\beta||_{2}^{2}+\alpha||\beta||_{1}\right],$

where $\lambda>0$ and $0\leq\alpha\leq 1$ are tuning parameters. Here
$\ell(a,b)$ is the negative log-likelihood for an observation e.g. for the
linear Gaussian case it is $\frac{1}{2}(a-b)^{2}$ and for logistic regression
it is $-a\cdot b+\log(1+e^{b})$.

The elastic net has been shown to be especially useful when the variables are
dependent [Zou and Hastie, 2005] (LD structure), which is the case of genetic
data and especially with bacterial genetic data. It is controlled by $\alpha$
that balances between lasso ($\alpha=1$) and ridge ($\alpha=0$). The tuning
parameter $\lambda$ controls the overall strength of the penalty and we use
10-fold cross-validation to choose suitable value for $\lambda$. Elastic net
is implemented in the software ’pyseer’ [Lees et al., 2018] focusing on
bacterial GWAS data and has been shown to outperform some state-of-the-art
methods in bacterial GWAS [Saber and Shapiro, 2020, Lees et al., 2020].

To obtain a 95% confidence intervals, we propose to use a post-selection
strategy by fitting the OLS estimator on the selected covariates and obtain
the 95% confidence intervals for the estimated effect sizes. This is followed
from [Belloni et al., 2013].

#### Heritability estimation

Heritability estimates are performed using the Boosting Heritability approach
in [Mai et al., 2021] where heritability estimation of commonly used methods
is improved by using a sample splitting strategy.

##  Appendix B Supplementary additional results

Figure 5:  Sample Correlation matrix of phenotypes in Maela data using binary
(left) and zone size (right) data

Figure 6: Histogram of zone size phenotypes in Maela data Table 1: Binary
Maela data. Heritability in the diagonal, genetic corelation in lower
diagonals, genetic related in upper diagonals

Chloramphenicol Erythromycin Penicillin Tetracycline Co-trimoxazole
carriage.duration Chloramphenicol 0.763 0 0 0 0 -0.009 (0.680, 0.838) (0, 0)
(0, 0) (0, 0) (0, 0) (-0.091, 0.011) Erythromycin -0.017 0.915 0.004 0 0 0
(-0.185, 0.707) (0.889, 0.941) (-0.002, 0.009) (0, 0) (-0.009, 0.001) (-0.009,
0.001) Penicillin -0.108 0.586 0.829 0 -0.671 0 (-0.326, 0.458) (-0.999,
0.545) (0.800, 0.851) (0, 0) (-0.999, 0.210) (0, 0) Tetracycline 0.025 0.441
0.261 0.891 0.007 0 (-0.156, 0.191) (-0.999, 0.543) (-0.383, 0.371) (0.873,
0.908) (0.001, 0.031) (0, 0) Co-trimoxazole -0.127 0.448 0.571 0.284 0.638 0
(-0.343, 0.370) (-0.999, 0.510) (-0.999, 0.934) (-0.430, 0.371) (0.565, 0.709)
(0, 0) carriage.duration 0.029 0.436 0.241 0.253 0.299 0.459 (-0.131, 0.145)
(-0.002, 0.882) (-0.091, 0.783) (-0.079, 0.257) (-0.080, 0.733) (0.356, 0.547)

Table 2: Binary Maela data. genetic covariance in the lower diagonals, co-
variance of causal effects in upper diagonals.

Chloramphenicol Erythromycin Penicillin Tetracycline Co-trimoxazole
carriage.duration Chloramphenicol 0 0 0 0 -0.080 (0, 0) (0, 0) (0, 0) (0, 0)
(-0.541, 0.264) Erythromycin -0.001 0.003 0 -0.001 0 (-0.126, 0.115) ()
(-0.001, 0.015) (0, 0) (-0.005, 0.001) (0, 0) Penicillin -0.010 0.095 0 -1.501
0 (-0.170, 0.144) (-0.480, 0.831) () (0, 0) (-4.433, 1.021) (0, 0)
Tetracycline 0.002 0.068 0.045 0.015 0 (-0.122, 0.124) (-0.228, 0.477)
(-0.545, 0.923) (0.005, 0.029) (0, 0) Co-trimoxazole -0.009 0.055 0.078 0.037
0 (-0.816, 0.671) (-0.686, 0.890) (-4.227, 5.074) (-0.782, 1.041) (0, 0)
carriage.duration 0.006 0.154 0.095 0.096 0.090 (-0.312, 0.364) (-0.001,
0.529) (-0.508, 0.783) (-0.203, 0.581) (-0.563, 0.733)

Table 3: Zone sizes Maela data. (unnormalized) genetic covariance in the lower
diagonals, (unnormalized) covariance of causal effects in upper diagonals.

Chloramphenicol Erythromycin Penicillin Tetracycline Co-trimoxazole
carriage.duration Chloramphenicol 0.447 0.166 0.338 0.219 0 (0.269, 0.657)
(0.099, 0.249) (0.207, 0.498) (0.163, 0.282) (0, 0) Erythromycin .026 0.437
0.848 0.383 0 (-.0298, .0716) (0.276, 0.634) (0.589, 1.161) (-0.228, 1.205)
(0, 0) Penicillin -0.030 0.542 0.010 1.040 0 (-0.109, 0.037) (0.225, 0.978)
(-4.418, 4.194) (0.536, 1.756) (0, 0) Tetracycline 0.029 0.434 0.096 0.482 0
(0.019, 0.042) (0.368, 0.507) (0.056, 0.146) (0.324, 0.669) (0, 0) Co-
trimoxazole -0.008 0.728 1.056 0.501 0.031 (-0.378, 0.248) (-1.375, 3.515)
(-3.503, 6.090) (-0.281, 1.499) (-0.065, 0.142) carriage.duration 0.017 0.391
0.290 0.187 0.423 (-0.050, 0.055) ( 0.200, 0.807) (-1.341, 0.619) ( 0.184,
0.651) (-2.356, 1.056)

Figure 7: Left: sample correlation matrix of phenotypes in Streptococcus
pneumoniae invasiveness data. Right: sample correlation matrix of log(MIC)
phenotypes in Neisseria gonorrhoeae data Table 4: Streptococcus pneumoniae
data. genetic covariance in the lower diagonals, covariance of causal effects
in upper diagonals.

Chloramphenicol Erythromycin Penicillin Tetracycline Co-trimoxazole
invasiveness Chloramphenicol 0.110 0.063 0.296 0.240 0.172 (0.072, 0.154)
(-0.097, 0.267) (0.145, 0.485) (0.127, 0.374) (0.119, 0.272) Erythromycin
0.011 0.137 0.032 0.093 0.082 (0.003, 0.027) (-0.043, 0.123) (0.061, 0.238)
(0.026, 0.187) (0.037, 0.142) Penicillin 0.002 0.007 0.096 0.853 0.126
(-0.023, 0.027) (-0.006, 0.029) (-5.122, 4.327) (-10.192, 9.067) (-0.356,
0.733) Tetracycline 0.003 0.002 0.001 0.083 0.018 (-0.003, 0.009) (-0.002,
0.008) (-0.293, 0.278) (-0.865, 1.068) (-1.013, 0.982) Co-trimoxazole 0.004
0.005 0.010 0.003 0.177 (-0.023, 0.028) (-0.013, 0.020) (-0.712, 0.714)
(-0.024, 0.030) (0.022, 0.453) invasiveness 0.001 0.001 0.001 0.002 0.002
(-0.027, 0.031) (-0.016, 0.019) (-0.023, 0.033) (-0.029, 0.031) (-0.002,
0.009)

Table 5: Ecoli data. genetic covariance in the lower diagonals, covariance of
causal effects in upper diagonals.

Amoxicillin Cefotaxime Ceftazidime Cefuroxime Ciprofloxacin Gentamicin
Amoxicillin 0 0 0 -0.0002 0 (0, 0) (0, 0) (0, 0) (-0.007, 0.005) (0, 0)
Cefotaxime -0.010 0.046 0.032 -0.013 0.007 (-0.297, 0.246) (-0.119, 0.244)
(-0.093, 0.156) (-0.073, 0.040) (-0.031, 0.038) Ceftazidime 0.0 0.027 0.005
0.010 0.006 (-0.181, 0.181) (-2.693, 2.991) (-0.046, 0.091) (-0.037, 0.060)
(-0.016, 0.037) Cefuroxime -0.009 0.041 0.022 0.044 0.006 (-0.179, 0.131)
(-1.969, 2.339) (-1.111, 1.268) (-0.059, 0.126) (-0.015, 0.028) Ciprofloxacin
-0.008 0.045 0.029 0.055 0.100 (-0.032, 0.006) (-0.396, 0.493) (-0.205, 0.268)
(-0.195, 0.298) (0.046, 0.159) Gentamicin 0.0005 0.017 0.015 0.019 0.037
(-0.235, 0.239) (-3.280, 3.530) (-1.922, 2.035) (-1.997, 2.143) (-0.390,
0.474)

Table 6:  Neisseria gonorrhoeae data. genetic covariance in the lower diagonals, covariance of causal effects in upper diagonals. | Azithromycin | Cefixime | Ciprofloxacin | Penicillin | Tetracycline
---|---|---|---|---|---
Azithromycin |  | 0 | 0 | -0.040 | 0
| (0, 0) | (0, 0) | (-1.834, 1.641) | (0, 0)
Cefixime | -0.019 |  | -0.507 | 0.106 | 0
(-0.175, 0.184) |  | (-2.413, 1.125) | (-0.059, 0.487) | (0, 0)
Ciprofloxacin | -0.865 | -0.173 |  | 0.047 | 0.301
(-5.776, 2.457) | (-0.904, 0.341) |  | (-0.194, 0.308) | (0.162, 0.480)
Penicillin | -0.009 | 0.028 | 2.259 |  | -0.037
(-6.299, 5.817) | (-0.408, 0.351) | (-2.824, 10.803) |  | (-0.466, 0.402)
Tetracycline | -0.153 | 0.004 | 1.469 | 0.482 |
(-0.725, 0.146) | (-0.187, 0.243) | (0.163, 2.770) | (0.205, 0.832) |


# Chemoreception and chemotaxis of a three-sphere swimmer

Stevens Paz  Departamento de Matemáticas, Universidad del Valle, Calle 13 No.
100-00, Cali, Colombia.  Roberto F. Ausas  Instituto de Ciências Matemáticas e
de Computação, Universidade de São Paulo, São Carlos, Brazil.  Juan P.
Carbajal  Ostschweizer Fachhochschule, Institut für Energietechnik IET,
Rapperswil, Switzerland.  Gustavo C. Buscaglia  Instituto de Ciências
Matemáticas e de Computação, Universidade de São Paulo, São Carlos, Brazil.

###### Abstract

The coupled problem of hydrodynamics and solute transport for the Najafi-
Golestanian three-sphere swimmer is studied, with the Reynolds number set to
zero and Péclet numbers (Pe) ranging from 0.06 to 60. The adopted method is
the numerical simulation of the problem with a finite element code based upon
the FEniCS library.

For the swimmer executing the optimal locomotion gait, we report the Sherwood
number as a function of Pe in homogeneous fluids and confirm that little gain
in solute flux is achieved by swimming unless Pe is significantly larger than
10.

We also consider the swimmer as an learning agent moving inside a fluid that
has a concentration gradient. The outcomes of Q-learning processes show that
learning locomotion (with the displacement as reward) is significantly easier
than learning chemotaxis (with the increase of solute flux as reward). The
chemotaxis problem, even at low Pe, has a varying environment that renders
learning more difficult. Further, the learning difficulty increases severely
with the Péclet number. The results demonstrate the challenges that natural
and artificial swimmers need to overcome to migrate efficiently when exposed
to chemical inhomogeneities.

##  1 Introduction

Locomotion at microscopic scales is a fundamental process in several areas of
biology, such as the study of the microbial life cycle, planktonic ecology,
algal blooms, cell migration in tissue, fertilization and embryo development,
and immune response. Liquid environments are the norm for microscopic
creatures, which swim and crawl in many different and ingenious ways. Interest
in the Fluid Mechanics of cell and micro-organism motion motivated early work
by Taylor [1], Lighthill [2] and Purcell [3], among others. Over the years,
significant advances have been made that can be found in recent comprehensive
reviews, such as those by Guasto et al. [4], Elgeti et al. [5], Lauga [6] and
Goldstein [7, 8].

In the last decade, interest in microscopic locomotion has expanded beyond
natural processes toward artificial ones. In fact, the technology to deploy
self-propelled micro-robots is becoming increasingly available or, at least,
achievable [9, 10, 11, 12, 13, 14].

The swimming strategies employed by natural and synthetic micro-organisms are
highly diverse. To gain understanding of the basic mechanisms, theoretical
studies have frequently relied on simplified models. Perhaps the simplest
model is the squirmer, which consists of a spheroidal particle with a given
slip velocity over its surface and was introduced by Lighthill [15] and Blake
[16] (see also Stone and Samuel [17]). Another simple early model is the
Taylor sheet [1] (see also Childress [18]), a flexible sheet that swims by
deforming transversely according to a moving wave.

In this study we address the three-sphere swimmer model, introduced by Najafi
and Golestanian [19, 20] inspired by an earlier three-link mechanism proposed
by Purcell [3]. This swimmer consists of three aligned spheres, with two
contractible arms (of zero thickness) joining them (see Fig. 1). It exhibits
quite large surface motions because of the variable arm length with the
simplicity of one-dimensional kinematics and spherical constituents. Because
of its mathematical simplicity, the three-sphere swimmer has become a popular
model [21, 22, 23, 24, 25, 26] with the additional advantage of having been
realized experimentally [27].

Our focus is on the diffusive transport of dissolved species to (or from) the
surface of the swimmer. This process is known as chemoreception. When the
dissolved species is a nutrient, chemoreception governs the nutrient uptake of
the organism and is central to its biological viability. The first studies
[28, 29] considered a sphere in steady motion. Later, a spherical squirmer
(self-propelled) was analyzed by Magar et al. [30] and Magar and Pedley [31]
and a suspension of spherical squirmers was simulated by Ishikawa et al. [32].

The chemoreception of kinematically more complex swimmers, such as the three-
sphere swimmer model, has not been reported in the literature. With the
increased interest in three-sphere swimmers and artificial swimmers in
general, the first goal of this contribution is to make available a first
study of the chemoreception of these swimmers for several swimming gaits and
Péclet numbers. We adopt a numerical approach based on the finite element
method which is an evolution of published work [33, 34]. Previous approaches
have mainly relied on Oseen’s approximation of the velocity field, which
requires the three spheres to be very far apart. One of the few published
similar works considers a biflagellated sphere [35], with the flagella moving
according to the swimming gates of the alga Chlamydomonas.

Our second goal concerns chemotaxis, that is, the ability of the swimmer to
move towards regions of higher solute (nutrient) concentration based on
chemoreception data. Recent attempts have shown that artificial swimmers can
learn locomotion by reinforcement learning algorithms [14, 36, 26, 37, 38,
39]. Learning chemotaxis was considered by Hartl et al. [40], who applied a
genetic learning algorithm to the three-sphere swimmer and showed that it can
successfully navigate static and time-dependent chemical environments.
However, they considered a simplified fluid dynamical model that did not
account for the swimmer’s impact on the concentration field. In this
contribution we report on some learning experiments with realistic interaction
of the swimmer with the environment. A basic Q-learning algorithm [41] is used
to process data from random sequences of actions. Then, from looking at the
learning outcomes, the increased difficulty of learning chemotaxis over
learning locomotion becomes apparent. The difficulty grows severely for Péclet
numbers above 1.

##  2 Governing equations

We consider the one-dimensional Najafi-Golestanian swimmer [19], which
consists of three spheres of radius $R$ joined by two extensible virtual (zero
thickness) arms (see Fig. 1). The variable lengths $L_{1}(t)$ and $L_{2}(t)$
of the arms, which are given functions of time, produce net displacements (in
the laboratory frame) $X_{1}(t)$, $X_{2}(t)$ and $X_{3}(t)$ of the centers of
the spheres and thus of the center of mass of the swimmer, $X(t)$. Notice
that, since $X=(X_{1}+X_{2}+X_{3})/3$,

$X_{1}=X-\frac{2L_{1}}{3}-\frac{L_{2}}{3},\quad
X_{2}=X+\frac{L_{1}}{3}-\frac{L_{2}}{3},\quad
X_{3}=X+\frac{L_{1}}{3}+\frac{2L_{2}}{3}~{}.$

The centers $\mathbf{C}_{i}$ ($i=1,2,3$) of the spheres are assumed lying
along the $x_{1}$ axis, i.e., $\mathbf{C}_{i}=(X_{i},0,0)^{T}$. The domain
occupied by the solid spheres is thus

$\mathcal{B}(X,L_{1},L_{2})=\mathcal{B}_{1}\cup\mathcal{B}_{2}\cup\mathcal{B}_{3}~{},\qquad\text{with
}\quad\mathcal{B}_{i}=\left\\{\mathbf{x}\in\mathbb{R}^{3}~{}|~{}\|\mathbf{x}-\mathbf{C}_{i}\|\leq
R\right\\}$

The swimmer evolves inside a large domain $\Omega$ with boundary conditions
mimicking an infinite medium. The fluid domain is
$\Omega_{f}\left(X,L_{1},L_{2}\right)=\Omega\setminus\mathcal{B}\left(X,L_{1},L_{2}\right)$,
filled with an inertialess incompressible Newtonian fluid of viscosity $\mu$.
The Cauchy stress of the fluid is thus given by

$\boldsymbol{\sigma}=-p\,\mathbf{I}+\mu\,\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right)~{},$
(1)

where $p$ and $\mathbf{u}$ are the fluid’s pressure and velocity fields,
respectively, and $\mathbf{I}$ is the identity matrix.

Figure 1: Three-spheres swimmer considered by Najafi and Golestanian [19].

Assuming no external forces acting on the swimmer, the hydrodynamic problem
reads:

Given $X(0)$, $L_{1}(t)$ and $L_{2}(t)$, for $t>0$ find $X(t)$,
$\mathbf{u}(\mathbf{x},t)$ and $p(\mathbf{x},t)$ such that

$\displaystyle\nabla\cdot\boldsymbol{\sigma}$ $\displaystyle=$ $\displaystyle
0~{},\qquad\text{in }\Omega_{f}(t)$ (2) $\displaystyle\nabla\cdot\mathbf{u}$
$\displaystyle=$ $\displaystyle 0~{},\qquad\text{in }\Omega_{f}(t)$ (3)
$\displaystyle\int_{\partial\mathcal{B}(t)}\boldsymbol{\sigma}\cdot\mathbf{\check{n}}~{}dS$
$\displaystyle=$ $\displaystyle 0~{},$ (4) $\displaystyle\mathbf{u}$
$\displaystyle=$
$\displaystyle\left(\frac{dX}{dt}-\frac{2}{3}\frac{dL_{1}}{dt}-\frac{1}{3}\frac{dL_{2}}{dt}\right)\hat{\mathbf{e}}~{},\qquad\text{on
}\partial\mathcal{B}_{1}$ (5) $\displaystyle\mathbf{u}$ $\displaystyle=$
$\displaystyle\left(\frac{dX}{dt}+\frac{1}{3}\frac{dL_{1}}{dt}-\frac{1}{3}\frac{dL_{2}}{dt}\right)\hat{\mathbf{e}}~{},\qquad\text{on
}\partial\mathcal{B}_{2}$ (6) $\displaystyle\mathbf{u}$ $\displaystyle=$
$\displaystyle\left(\frac{dX}{dt}+\frac{1}{3}\frac{dL_{1}}{dt}+\frac{2}{3}\frac{dL_{2}}{dt}\right)\hat{\mathbf{e}}~{},\qquad\text{on
}\partial\mathcal{B}_{3}$ (7)

Equations (2)-(3) are the well known Stokes equations, (4) is the force-free
condition where $\mathbf{\check{n}}$ is the unit normal outwards of
$\Omega_{f}$, and (5)-(7) impose the no-slip condition at the surface of each
sphere, being $\hat{\mathbf{e}}$ the unit positive direction of motion.

Coupled with the hydrodynamic problem is the solute transport problem, in
which the unknown is the concentration field $C(\mathbf{x},t)$. It reads:

Given the initial concentration field $C(\mathbf{x},t=0)$, for $t>0$ find
$C(\mathbf{x},t)$ satisfying

$\displaystyle\frac{\partial C}{\partial t}+\mathbf{u}\cdot\nabla
C-D\,\nabla^{2}C$ $\displaystyle=$ $\displaystyle 0~{},\qquad\qquad\text{in
}\Omega_{f}~{},$ (8) $\displaystyle D\nabla C\cdot\mathbf{\check{n}}$
$\displaystyle=$ $\displaystyle-\phi_{i}\,C~{},\qquad\text{on
}\partial\mathcal{B}_{i}~{}.$ (9)

with boundary conditions $C=C_{\infty}$ far away from the swimmer.

In (8), $D$ is the diffusivity of the solute. The so-called Robin boundary
condition (9) contains the parameters which govern the permeability of the
spheres’ surfaces to the solute. As $\phi_{i}\to\infty$, the Robin condition
tends to the Dirichlet condition $C=0$. Of crucial importance in our analysis
is the total solute flux $J$ to the swimmer, defined as

$J(t)=\int_{\partial\mathcal{B}(t)}D\nabla C\cdot\mathbf{\check{n}}~{}dS~{}.$
(10)

##  3 Numerical method: Description and validation

The numerical method is an adaptation of the Arbitrary Lagrangian-Eulerian
(ALE) finite element methods discussed in previous work by the authors in the
context of squirmers [33] and of flexible swimmers [34]. The current
implementation relies heavily on the FEniCS package [42].

The spatial interpolation of velocity and pressure adopts $P_{1}$ conforming
finite elements, stabilized by the well-known ASGS (or GLS) technique [43, 44,
45]. The interpolation of the concentration field is carried out with $P_{2}$
finite elements. Stabilization of advection is also performed by the ASGS
method. Large deformations of the mesh are mitigated by elastic mesh updates
plus frequent remeshing based on a mesh quality criterion.

The time stepping strategy follows a second-order predictor-corrector method
for the hydrodynamic problem and a second-order trapezoidal scheme for the
solute transport. More details can be found in Appendix A.

Several tests were performed to validate the numerical code:

a)

Displacement achieved by the swimming gait of the Najafi-Golestanian swimmer.
By modifying the length of one arm at a time, between the maximum length $W$
and the minimum one $W-w$, the swimmer advances a distance $\Delta X$ per
gait. Previous asymptotic and numerical estimates [19, 46] are compared in
Fig. 2 to our numerically obtained values, showing excellent agreement.

b)

Chemoreception of a sphere towed at constant velocity in a fluid. It is known
that for a towed sphere the solute flux $J$, for Dirichlet conditions
($\phi\to+\infty$), is well approximated by the expression proposed by Clift
et al. [47] and Karp-Boss et al. [48]

$\text{Sh}=\frac{J}{J_{0}}=\frac{1}{2}\left(1+\left(1+2\text{Pe}\right)^{\frac{1}{3}}\right)~{},$
(11)

where $J_{0}$ is the solute flux when towing velocity $V$ is zero (purely
diffusive transport) and Pe $=VR/D$ is the Péclet number. The ratio $J/J_{0}$
is known as Sherwood number Sh. Our numerical results were verified to agree
within 3% with Clift’s formula in the whole range $10^{-2}\leq\text{Pe}\leq
10^{2}$.

Figure 2: Asymptotic behaviour of normalized displacement $\Delta X/R$ when
$w/R\to 0$ for $3$-linked spheres swimmer. The $P1/P1$-stabilized solution
(asterisk marker) is compared against the formulae from Earl et al. [46] and
Najafi and Golestanian [19].

##  4 Solute flux toward a three-sphere swimmer

The baseline solute flux $J_{0}$ of a three-sphere swimmer is defined here as
the steady diffusive flux at its surface when the arms of the swimmer are
totally extended. In what follows we concentrate on the Dirichlet condition
case, i.e., all $\phi_{i}\to+\infty$. If $W\gg R$ the analytic solution
$C(r)=C_{\infty}(1-\frac{R}{r})$ for a sphere yields the baseline flux
$\hat{J}_{0}=12\pi RDC_{\infty}$. For smaller values of $W/R$ one gets smaller
values of $J_{0}$, as shown in Fig. 3.

Figure 3: $J_{0}/\hat{J}_{0}$ as a function of $W/R$.

In this study we consider the basic non-reciprocal swimming gait discussed by
Najafi and Golestanian [19], by which the arms move one at a time, taking
turns (i.e., starting from the extended position, contract arm 1, contract arm
2, extend arm 1, extend arm 2, repeat) as shown in Fig. 4 and Fig. 5. An
equivalent gait is obtained by reversing the sequence in time, the only
difference being that the direction of the movement is to the left.

Figure 4: Swimming gait for the three-spheres swimmer considered by Najafi and
Golestanian [19]. The sequence is from top to bottom and follows the
numeration introduced in Sec. 5.1 for states ($s$) and actions ($a$).

Figure 5: Swimming gait for the three-spheres swimmer considered by Najafi and
Golestanian [19]. Left: velocity magnitude and streamlines. Right: pressure
field and contour lines. In both cases, colors going from red to blue
indicating maximum to minimum quantity value.

When it’s time for one specific arm to move, it goes to complete extension
($L_{i}=W$) to complete contraction ($L_{i}=W-w$), if it was extended, or
viceversa, if it was contracted. Though the net displacement obtained by
executing the gait does not depend on the speed at which the sequence of
movements is made (because fluid inertia is negligible), the solute transport
is affected by the speed. The simplest case was adopted, in which all arm
extensions and contractions are performed at uniform, constant speed $S$. The
duration of the gait is thus $T=4w/S$, after which the swimmer recovers its
initial posture, but its position is affected by a translation $\Delta X$. The
computation of $\Delta X$ requires solving the hydrodynamic problem (2)-(7).
The values of $\Delta X$ obtained with our numerical code are compared to
reference ones published by Earl et al (2007) and to asymptotic approximations
in Fig. 2. The good agreement is an additional validation of the code used in
this work. The average velocity of the center of mass of the swimmer is
$V=S\Delta X/(4w)$. Let us denote $\Delta\mathbf{X}=(\Delta X,0,0)^{T}$. The
solution of the hydrodynamic problem is periodic in the swimmer’s frame and
thus only needs to be solved in time along one gait duration (say, between
$t=0$ and $t=T$), since

$\displaystyle X(t+T)$ $\displaystyle=$ $\displaystyle X(t)+\Delta X~{},$ (12)
$\displaystyle\mathbf{u}(\mathbf{x}+\Delta\mathbf{X},t+T)$ $\displaystyle=$
$\displaystyle\mathbf{u}(\mathbf{x},t)~{},$ (13) $\displaystyle
p(\mathbf{x}+\Delta\mathbf{X},t+T)$ $\displaystyle=$ $\displaystyle
p(\mathbf{x},t)~{}.$ (14)

The transport problem (8)-(9), on the other hand, depends upon the initial
condition and evolves in time as the immediate surroundings of the swimmer are
gradually depleted and the boundary layers produced by the swimmer’s movement
develop. Eventually, after some development time $\mathcal{T}$, the
concentration field also reaches a periodic evolution in the swimmer’s frame,
namely, for $t>\mathcal{T}$,

$C(\mathbf{x}+\Delta\mathbf{X},t+T)=C(\mathbf{x},t)~{}.$ (15)

This implies that $J(t)$ is also periodic, $J(t+T)=J(t)$. Our aim is to
provide data about the average solute flux

$\overline{J}=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}J(t)~{}dt~{}$ (16)

once periodic conditions are attained, i.e., $t_{0}>\mathcal{T}$. The average
solute flux when the swimmer is executing the swimming gait is denoted by
$J_{s}(\text{Pe})$.

Remark: The system necessarily arrives at a periodic state because the
hydrodynamic problem is totally decoupled from the transport problem, the
latter being linear and parabolic (albeit in a periodically-varying domain).
Nonlinearities by which the concentration field affects the hydrodynamics may
lead to more complex behaviors.

The parameters of the problem are $R$, $W$, $w$, $D$ and $S$. The far field
concentration $C_{\infty}$ is irrelevant, taken as unity. The relevant non-
dimensional variables are $W/R$, $w/R$ and the Péclet number

$\text{Pe}=\frac{SR}{D}~{}.$ (17)

We have considered $W/R=10$ combined with different values of $w/R=2,4$ and
$6$. Notice that in the last case the spheres are just one diameter apart when
the arm is in the contracted position.

Figure 6: Sherwood number as a function of Pe for $w/R=2,4$ and $6$.

The solution of the full problem (hydrodynamics+transport) yields the average
Sherwood number $J_{s}/J_{0}$ as a function of Pe for each case, as shown in
Fig. 6. It is worth pointing out that for very low Pe the solute flux is
maximal when the swimmer is immobile with both arms extended. The swimming
gait takes the arms out of the optimal position and thus results in
$J_{s}<J_{0}$, as evinced by Fig. 6 for Pe $<5$.

In Fig. 7 we present the transient fluxes $J(t)/J0$ starting from an initial
condition that is the steady-state solution corresponding to the initial
posture of the swimmer. We include, in Fig. 8, some sample snapshots of the
concentration wake of the swimmer at different times.

Figure 7: Left: Time histories of $J(t)$ up to non-dimensional time $t=100$
after starting the swimming gait at different Péclet numbers. Right: Long term
behavior of the flux averaged over the swimming gait (four actions) up to time
$t=500$.

$tS/w=0$

$tS/w=240$

$tS/w=500$

Figure 8: Snapshots of the concentration wake for different times and Péclet
numbers (contour lines in black color and colors going from red to blue
indicating maximum to minimum concentration). Top: initial condition for all
cases. Center: time $tS/w=240$ and Pe $=60,6,0.6$ (top to bottom,
respectively). Bottom: time $tS/w=500$ and Pe $=60,6,0.6$ (top to bottom,
respectively)

##  5 Q-Learning chemotaxis

###  5.1 The reinforcement learning framework

One can view the swimmer as an agent in discrete time performing a decision
process [49]. The time unit is $w/S$, the time it takes to contract an
extended arm or vice versa. In its simplest formulation, the agent has four
possible states: both arms contracted (state 1), arm 1 contracted and arm 2
extended (state 2), arm 1 extended and arm 2 contracted (state 3) and both
arms extended (state 4). The sequence of states that corresponds to the
swimming gait (toward the right, i.e. $x_{1}>0$) is 1-3-4-2-1-3-$\ldots$. The
immobile extended mode is a sequence 4-4-$\ldots$.

We define that the agent must move exactly one arm in each time interval, from
the fully extended position to the fully contracted position or vice versa.
The only decision to be taken at each discrete instant $t$ is whether to move
arm 1 (action 1) or arm 2 (action 2). At each of the four possible states
there are just two possible actions. This kinematic restriction on the agent
rules out the possibility of the swimmer remaining still, in any position. The
swimming gait to the right, starting at state 4, corresponds to the sequence
of actions 1-2-1-2-$\ldots$ .

The agent defined above can be made to follow a policy, i.e., a rule that to
each state $s_{t}\in\\{1,2,3,4\\}$ associates a probability
$\pi_{t}(s_{t},a=1)$ of taking action $a=1$ at time $t$ (and thus
$\pi_{t}(s_{t},a=2)=1-\pi_{t}(s_{t},a=1)$ of taking action $a=2$). We only
consider here time-independent policies, so that we drop the subscript $t$
from $\pi_{t}$ hereafter. A totally random policy selects, at each decision
time, action 1 or 2 randomly with 50% probability each (i.e.,
$\pi(s_{t},1)=\pi(s_{t},2)=0.5$). The swim-to-the-right gait is a totally
deterministic policy that takes action 1 if $s_{t}=1$ or $s_{t}=4$ and takes
action 2 if $s_{t}=2$ or $s_{t}=3$. This corresponds to $\pi(1,1)=\pi(4,1)=1$
and $\pi(2,1)=\pi(3,1)=0$. This particular policy is denoted by $\pi^{*}$.
Other policies are of course possible, and the identification of a convenient
policy is referred to as learning.

Whichever the policy, each action causes the agent to interact with the fluid.
Assume action $a_{t}$ is taken at time $t$, when the agent is at state $s_{t}$
in an environment described by the variable $\theta_{t}$. After the action,
the state of the agent is altered to another state $s_{t+1}$ and the
environment changes to $\theta_{t+1}$. In our setting, there is only one
$s_{t+1}$ possible as a consequence of each action, and it is independent of
the environment.

The environment variable $\theta_{t}$ represents all the relevant information
about the environment at time $t$. To define it more precisely, we consider
our agents moving inside a domain $\Omega$ that is fixed in the laboratory
frame. The geometry of the fluid environment (the domain $\Omega_{f}$
discussed in Section 2) is fully determined, at time $t$, by $X_{t}=X(t)$,
$L_{1}(t)$ and $L_{2}(t)$. The environment variable here is $X_{t}$, since
$L_{1}(t)$ and $L_{2}(t)$ are state variables. Looking now at the solute
transport problem, the only additional information that affects the evolution
of the system between $t$ and $t+1$ is the initial concentration field
$C_{t}=C(\mathbf{x},t)$ (with $\mathbf{x}$ arbitrary in $\Omega_{f}(t)$). The
environment variable $\theta_{t}$ is thus defined as
$\theta_{t}=(X_{t},C_{t})$. By solving the hydrodynamic+transport problem from
time $t$ to time $t+1$ one can update the environment variables. Formally, we
write this as

$\theta_{t+1}=\Theta(s_{t},a_{t},\theta_{t})~{},$

where the update function $\Theta$ is not explicitly known.

The learning problem we consider is that of finding a policy $\pi^{*}$ that
maximizes some additive total reward

$R(\pi)=\sum_{\tau>0}\gamma^{\tau}r_{\tau}~{},$ (18)

where $\gamma\in(0,1)$ is the discount factor and $r_{\tau+1}$ is the
instantaneous reward obtained by performing the action $a_{\tau}$ selected by
the policy $\pi$. As $\gamma$ increases toward 1, the total reward gives more
importance to the long term consequences of current actions.

###  5.2 Rewards associated to chemotaxis

Consider the 3-sphere agent immersed in a fluid where there exists a constant
gradient of solute concentration, i.e., far away from the agent the
concentration field is given by $C(x,y,z,t)=g\,x$, where $g>0$. Assume that
the agent wants to maximize its access to the solute, then it is easily seen
that its optimal policy is nothing but the swimming gait (to the right). In
terms of an additive reward this can be set in several ways, of which we
consider the following three:

R1)

Displacement reward: In this case the instantaneous reward is simply

$r_{t+1}=X_{t+1}-X_{t}~{},$ (19)

so that $\sum_{t=1}^{N}r_{t}=X(N)-X(0)$. Notice that this reward is purely
hydrodynamic (it is independent of the concentration field). It corresponds to
”learning to swim”, or ”learning locomotion”, rather than ”learning
chemotaxis”. Tsang et al. [14] have shown that the Q-learning algorithm [41]
successfully identifies optimal policies for R1. They used a simplified
hydrodynamic model (Oseen approximation) that is only valid when $W-w\gg R$.
Our hydrodynamic model (2)-(7), being exact, is not subject to this condition.
In fact, we show below that Q-learning is also very efficient at learning to
swim when $W-w=4R$.

The rewards we consider next are true chemotactic rewards and, unlike R1, they
are observable to the agent without knowing its relative position, w.r.t. the
lab. frame. Our aim is to evaluate the difficulty of learning chemotaxis as
compared to that of learning locomotion.

R2)

Intake accumulation: This reward is the total intake of solute during the
interval $[t,t+1]$,

$r_{t+1}=\int_{t}^{t+1}J(\tau)~{}d\tau~{},$ (20)

so that $\sum_{t=1}^{N}r_{t}=\int_{0}^{N}J(\tau)~{}d\tau$. It expresses the
drive to maximize the total intake along the evolution of the system.

R3)

Intake increase: This reward is defined as

$r_{t+1}=J(t+1)-J(t)~{},$ (21)

so that $\sum_{t=1}^{N}r_{t}=J(N)-J(0)$. In drives the agent towards
maximizing the solute intake rate at the end of the evolution.

The consideration of the solute transport variables in rewards R2 and R3
increases the complexity of the agent-environment problem, because
$C(\cdot,t)$ is a field and thus high-dimensional.

###  5.3 The Q-learning algorithm

Q-learning is based on the value-iteration algorithm, which provably (under
suitable assumptions) converges to a solution of the Bellman optimality
equation.

The algorithm successively updates an action-value matrix $Q(s,a)$ which
represents the maximum total reward that can be expected if action $a$ is
taken when the agent is at state $s$.

An experience is defined as a 4-tuple $e_{t+1}=(s_{t},a_{t},s_{t+1},r_{t+1})$
consisting of the agent’s state prior to the action, the action taken, the
agent’s state after the action, and the instantaneous reward collected. Given
a set of experiences $e_{1},e_{2},\ldots,e_{N}$, and initializing $Q_{0}$ with
zeroes, $Q_{t}$ is updated as

$Q_{t+1}(s_{t},a_{t})=(1-\alpha)Q_{t}(s_{t},a_{t})+\alpha\,\left(r_{t+1}+\gamma\,\max_{a\in\\{1,2\\}}Q_{t}(s_{t+1},a)\right)~{},$
(22)

where $\alpha$ is the so-called learning rate. After the $N$ updates, one
obtains matrix $Q_{N}$. The learned policy $\pi_{\text{\tiny QL}}$ is totally
greedy in $Q_{N}$, i.e.,

$\pi_{\text{\tiny
QL}}(s,a)=0\quad\forall(s,a)~{}\mbox{s.t.}~{}Q_{N}(s,a)<\max_{a^{\prime}}Q_{N}(s,a^{\prime})~{}.$

For each set of experiences, the learned policy $\pi_{\text{\tiny QL}}$
depends on the two parameters $\gamma$ and $\alpha$.

Q-learning is an off-policy algorithm, in the sense that the set of
experiences can be taken from any evolution of the system that repeatedly
visits all possible states. We generate the experience sets with random
policies, as will be described later on. Along the generation of each set we
store rewards R1, R2 and R3 for each experience. This allows us to compare the
learning outcomes of the Q-learning algorithm for all rewards and different
environments (Pe number) over the same experience set.

###  5.4 Numerical experiments

We consider an agent with geometrical parameters $W/R=10$ and $w/R=6$ that is
immobile at state $s_{0}$ and center-of-mass position $X_{0}$ for a long time
in a homogeneous environment. The concentration field has reached, at time
$t=0$, its steady state with value $C_{\infty}(x,y,z)=g\,x$ far from the agent
and 0 at its surface. At $t=0$ we start the decision process following a
totally random policy and solve the hydrodynamic+transport problem along $N$
actions until time $t=N$, storing the sequence of experiences corresponding to
each reward. The solution is computed numerically with the FEniCS-based
solver, using a time step $\delta t=0.1$. The spatial resolution ensures that
the mesh size close to the spheres is about $0.05R$.

In this way, long sequences of actions were computed and stored for random
motions at different Péclet numbers ranging from Pe $=0.06$, which is fully
diffusion dominated, to Pe $=60$. The instantaneous rewards $r_{t}$
corresponding to the first $50$ actions of some of the recorded sequences are
plotted in Fig. 9, together with the position $X_{t}$.

Figure 9: Rewards (blue) and normalized center of mass position (red) for $50$
actions for a realization with Pe $=0.06$. The red vertical scale is the same
for all cases.

Let us first discuss the results obtained at Pe $=0.06$. This is the case with
lowest complexity, since the concentration field does not depend on the
history and evolves quasi-statically. We divided the generated sequence of
6000 experiences into $12$ sets of $500$ consecutive experiences and ran the
Q-learning process on each experience set for different values of $\alpha$ and
$\gamma$. The resulting policy $\pi_{\text{\tiny QL}}$ was compared to
$\pi^{*}$, classifying the learning process as successful if $\pi_{\text{\tiny
QL}}=\pi^{*}$.

The heat maps in Fig. 10 show the success rates obtained for each reward. For
reward R1 (displacement), the success rate is $100\%$ for all tested values of
$\alpha$ and $\gamma$. For reward R2 (intake accumulation) the situation is
quite the opposite, since the success rate is zero for almost every pair
$(\gamma,\alpha)$. In the case of reward R3 (intake increase), values of
$\gamma\geq 0.9$ show high success rates, but any $\gamma\leq 0.7$ fails
consistently. This suggests that, even at low Pe, learning chemotaxis is more
difficult than learning locomotion. This can be explained by looking at the
top row of Figure 11, where the instantaneous rewards $r_{t+1}$ of actions
$a=1$ and $a=2$ are plotted against $X_{t}$ for the state $s_{t}=2$. Reward R1
does not depend on $X$, it does not vary along the evolution. The other
rewards are affine on $X_{t}$, which makes them evolve over time and increases
the difficulty of learning.

The sharp constancy of $\Delta X$ for each action is an additional
verification of the accuracy of the numerical model. Though the mesh is
different each time the action is executed, $\Delta X$ does not change.

Rewards R2 and R3 show some small noise on top of their affine dependence on
$X_{t}$. This is an effect of the non-zero Péclet number. If we now look at
the other rows of Figure 11, we see that the variance of the reward for each
action and each $X_{t}$ increases significantly with Pe, for both R2 and R3.
Already at Pe $=6$ the rewards of the two actions overlap for the same
$X_{t}$, making it even more difficult to predict the convenience of one
action over the other.

The consequences on the success rate of Q-learning at chemotaxis are as would
be expected. We extracted 10 batches of 500 consecutive experiences from the
the sequences of 1000 experiences collected at each Pe $\geq$ 0.6. In Fig. 12
we show the success rates over the ten batches as heat maps over the
parameters $\gamma$ and $\alpha$ as before.

The displacement reward is unaffected by Pe and thus remains at 100 % success
rate throughout the maps. Reward R2, that performs very poorly at Pe = 0.06,
improves somewhat for higher Pe, outperforming R3 at Pe = 6 and 10. Reward R3
exhibits high success rates at low Pe but they decline systematically as Pe
increases to 6. Some improvement is then observed at Pe = 60. These
observations are summarized in Fig. 13, which shows box plots of the total
success rate (summed over $\gamma$ and $\alpha$) as function of Pe,
aggregating the data corresponding to the ten learning batches of experiences.

Figure 10: Heat maps of percentage of learning success obtained from 12 sets
of 500 consecutive experiences for the data set with Pe = 0.06.

Pe $=0.06$

Pe $=0.6$

Pe $=6.0$

Pe $=60.0$

Figure 11: Instantaneous rewards plotted as function of the center-of-mass
position $X$ for different values of Pe. Only experiences that start at state
$s_{t}=2$ are included, with different colors for the cases $a_{t}=1$ (blue)
and $a_{t}=2$ (yellow). Figure 12: Heat maps of percentage of learning success
obtained from $10$ sets of $500$ consecutive experiences for the data sets
corresponding, from top to bottom, to Pe = 0.06, 0.6, 1, 6, 10 and 60. Figure
13: Box plot of the total success rate (summed over $\gamma$ and $\alpha$) for
each reward and different Péclet numbers. Each data point corresponds to one
of ten batchs of 500 experiences each.

##  6 Conclusions

We have examined the coupled problem of hydrodynamics and solute transport for
a three-sphere swimmer at Péclet numbers ranging from 0.06 to 60. The adopted
method is the numerical simulation of the problem with a finite element code
based upon the FEniCS library.

We first analyzed the swimming gait and produced previously unavailable data
on the dependence of the Sherwood number with the Péclet number in homogeneous
environment. They confirm that little gain in solute flux can be achieved by
swimming (with respect to staying still in the fully extended position) unless
Pe is significantly larger than 10.

We then turned to considering the swimmer as an agent moving inside a fluid
that has a concentration gradient (chemotaxis). Inside a positive gradient,
the swimmer/agent has potentially unbounded gain in solute flux if it swims
toward the right.

A Q-learning process was applied to sets of experiences collected along random
decisions of the swimmer/agent to evaluate the difficulty of learning
chemotaxis as compared to that of learning locomotion. The numerical evidence
shows that the latter is significantly easier than the former. The chemotaxis
problem even at low Pe has a varying environment that renders learning more
difficult. Further, the average success rate falls sharply for Pe$\geq 1$ and
the dispersion in the outcome increases.

The overall conclusion is that proprioception (information about $s_{t}$) and
chemoreception (information about $J(t)$) alone pose a learning problem for
chemotaxis that is quite hard. Several routes toward alleviating this
difficulty exist. Further sensory information can be added, such as the
fraction of the flux that is absorbed by each sphere. Information about the
concentration field away from the swimmer could possibly be made available by
antenna-like chemoreceptive sensors. Another route consists of considering
that the swimmer/agent has an ”innate” bias towards random motions different
from the totally random motions used in this study. This would induce a
preferential sampling of the space of sequences of actions that could improve
learning. These routes, set for future work, will help elucidate the
challenges that swimmers must overcome to learn chemotaxis.

##  7 Acknowledgements

The authors gratefully acknowledge the financial support received from the São
Paulo Research Foundation FAPESP/CEPID/CeMEAI grant 2013/07375-0, and from
Brazilian National Council for Scientific and Technological Development CNPq
grants 305599/2017-8 and 310990/2019-0. This research was carried out using
computational resources from the Cluster-Sala Jürgen Tischer, Department of
Mathematics, Universidad del Valle, Cali, Colombia. The authors are thankful
to the developers of the free and open source software used for this
publication.

## References

  * Taylor [1951] Geoffrey Ingram Taylor.  Analysis of the swimming of microscopic organisms.  _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ , 209(1099):447–461, 1951\.
  * Lighthill [1976] James Lighthill.  Flagellar hydrodynamics.  _SIAM review_ , 18(2):161–230, 1976.
  * Purcell [1977] Edward M Purcell.  Life at low reynolds number.  _American Journal of Physics_ , 45(1):3–11, 1977\.
  * Guasto et al. [2012] Jeffrey S Guasto, Roberto Rusconi, and Roman Stocker.  Fluid mechanics of planktonic microorganisms.  _Annual Review of Fluid Mechanics_ , 44:373–400, 2012.
  * Elgeti et al. [2015] Jens Elgeti, Roland G Winkler, and Gerhard Gompper.  Physics of microswimmers—single particle motion and collective behavior: a review.  _Reports on Progress in Physics_ , 78(5):056601, 2015.
  * Lauga [2016] Eric Lauga.  Bacterial hydrodynamics.  _Annual Review of Fluid Mechanics_ , 48:105–130, 2016.
  * Goldstein [2015] Raymond E Goldstein.  Green algae as model organisms for biological fluid dynamics.  _Annual Review of Fluid Mechanics_ , 47:343–375, 2015.
  * Goldstein [2016] Raymond E Goldstein.  Batchelor prize lecture fluid dynamics at the scale of the cell.  _Journal of Fluid Mechanics_ , 807:1–39, 2016.
  * Becker et al. [2003] Leif E Becker, Stephan A Koehler, and Howard A Stone.  On self-propulsion of micro-machines at low reynolds number: Purcell’s three-link swimmer.  _Journal of fluid mechanics_ , 490:15–35, 2003.
  * Peyer et al. [2012] Kathrin E Peyer, Arthur W Mahoney, Li Zhang, Jake J Abbott, and Bradley J Nelson.  Bacteria-inspired microrobots.  _Microbiorobotics_ , pages 165–199, 2012.
  * Gao et al. [2012] Wei Gao, Daniel Kagan, On Shun Pak, Corbin Clawson, Susana Campuzano, Erdembileg Chuluun-Erdene, Erik Shipton, Eric E Fullerton, Liangfang Zhang, Eric Lauga, et al.  Cargo-towing fuel-free magnetic nanoswimmers for targeted drug delivery.  _Small_ , 8(3):460–467, 2012.
  * Pak et al. [2015] On Shun Pak, Eric Lauga, C Duprat, and H Stone.  _Theoretical models of low-Reynolds-number locomotion_.  Royal Society of Chemistry, 2015.
  * Liao et al. [2019] Pan Liao, Liuxi Xing, Shiwu Zhang, and Dong Sun.  Magnetically driven undulatory microswimmers integrating multiple rigid segments.  _Small_ , 15(36):1901197, 2019.
  * Tsang et al. [2020] Alan Cheng Hou Tsang, Pun Wai Tong, Shreyes Nallan, and On Shun Pak.  Self-learning how to swim at low reynolds number.  _Physical Review Fluids_ , 5(7):074101, 2020.
  * Lighthill [1952] MJ Lighthill.  On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers.  _Communications on Pure and Applied Mathematics_ , 5(2):109–118, 1952.
  * Blake [1971] JR Blake.  A spherical envelope approach to ciliary propulsion.  _Journal of Fluid Mechanics_ , 46(1):199–208, 1971.
  * Stone and Samuel [1996] Howard A Stone and Aravinthan DT Samuel.  Propulsion of microorganisms by surface distortions.  _Physical Review Letters_ , 77(19):4102, 1996\.
  * Childress [1981] Stephen Childress.  _Mechanics of swimming and flying_ , volume 2.  Cambridge University Press, 1981.
  * Najafi and Golestanian [2004] Ali Najafi and Ramin Golestanian.  Simple swimmer at low reynolds number: Three linked spheres.  _Physical Review E_ , 69(6):062901, 2004.
  * Najafi and Golestanian [2005] Ali Najafi and Ramin Golestanian.  Propulsion at low reynolds number.  _Journal of Physics: Condensed Matter_ , 17(14):S1203, 2005.
  * Golestanian [2008] Ramin Golestanian.  Three-sphere low-reynolds-number swimmer with a cargo container.  _The European Physical Journal E_ , 25(1):1–4, 2008.
  * Golestanian and Ajdari [2008] Ramin Golestanian and Armand Ajdari.  Analytic results for the three-sphere swimmer at low reynolds number.  _Physical Review E_ , 77(3):036308, 2008.
  * Alouges et al. [2009] F Alouges, Antonio DeSimone, and Aline Lefebvre.  Optimal strokes for axisymmetric microswimmers.  _The European Physical Journal E_ , 28(3):279–284, 2009.
  * Pickl [2016] Kristina Pickl.  _Fully Resolved Simulation of Self-propelled Microorganisms on Supercomputers_.  PhD thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), 2016.
  * Nasouri et al. [2019] Babak Nasouri, Andrej Vilfan, and Ramin Golestanian.  Efficiency limits of the three-sphere swimmer.  _Physical Review Fluids_ , 4(7):073101, 2019.
  * Berti [2021] Luca Berti.  _Numerical methods and optimisation for micro-swimming_.  PhD thesis, Université de Strasbourg, 2021.
  * Grosjean et al. [2016] Galien Grosjean, Maxime Hubert, Guillaume Lagubeau, and Nicolas Vandewalle.  Realization of the najafi-golestanian microswimmer.  _Physical Review E_ , 94(2):021101, 2016.
  * Berg and Purcell [1977] Howard C Berg and Edward M Purcell.  Physics of chemoreception.  _Biophysical Journal_ , 20(2):193–219, 1977.
  * Wiegel [1983] Frederik W Wiegel.  Diffusion and the physics of chemoreception.  _Physics Reports_ , 95(5):283–319, 1983.
  * Magar et al. [2003] Vanesa Magar, Tomonobu Goto, and Timothy J Pedley.  Nutrient uptake by a self-propelled steady squirmer.  _The Quarterly Journal of Mechanics and Applied Mathematics_ , 56(1):65–91, 2003.
  * Magar and Pedley [2005] Vanesa Magar and Timothy J Pedley.  Average nutrient uptake by a self-propelled unsteady squirmer.  _Journal of fluid mechanics_ , 539:93–112, 2005.
  * Ishikawa et al. [2016] Takuji Ishikawa, Shunsuke Kajiki, Yohsuke Imai, and Toshihiro Omori.  Nutrient uptake in a suspension of squirmers.  _Journal of Fluid Mechanics_ , 789:481–499, 2016.
  * Paz and Buscaglia [2020] Stevens Paz and Gustavo C Buscaglia.  Simulating squirmers with volumetric solvers.  _Journal of the Brazilian Society of Mechanical Sciences and Engineering_ , 42(10):1–27, 2020.
  * Ausas et al. [2022] Roberto Federico Ausas, Cristian Guillermo Gebhardt, and Gustavo Carlos Buscaglia.  A finite element method for simulating soft active non-shearable rods immersed in generalized newtonian fluids.  _Communications in Nonlinear Science and Numerical Simulation_ , 108:106213, 2022.
  * Tam and Hosoi [2011] Daniel Tam and AE Hosoi.  Optimal feeding and swimming gaits of biflagellated organisms.  _Proceedings of the National Academy of Sciences_ , 108(3):1001–1006, 2011.
  * Misra et al. [2020] Dipendra Misra, Mikael Henaff, Akshay Krishnamurthy, and John Langford.  Kinematic state abstraction and provably efficient rich-observation reinforcement learning.  In _International conference on machine learning_ , pages 6961–6971. PMLR, 2020.
  * Alageshan et al. [2020] Jaya Kumar Alageshan, Akhilesh Kumar Verma, Jérémie Bec, and Rahul Pandit.  Machine learning strategies for path-planning microswimmers in turbulent flows.  _Physical Review E_ , 101(4):043110, 2020.
  * Qiu et al. [2022] Jingran Qiu, Navid Mousavi, Kristian Gustavsson, Chunxiao Xu, Bernhard Mehlig, and Lihao Zhao.  Navigation of micro-swimmers in steady flow: The importance of symmetries.  _Journal of Fluid Mechanics_ , 932, 2022.
  * Muiños-Landin et al. [2021] Santiago Muiños-Landin, A Fischer, V Holubec, and F Cichos.  Reinforcement learning with artificial microswimmers.  _Science Robotics_ , 6(52):eabd9285, 2021.
  * Hartl et al. [2021] Benedikt Hartl, Maximilian Hübl, Gerhard Kahl, and Andreas Zöttl.  Microswimmers learning chemotaxis with genetic algorithms.  _Proceedings of the National Academy of Sciences_ , 118(19), 2021.
  * Watkins and Dayan [1992] Christopher JCH Watkins and Peter Dayan.  Q-learning.  _Machine Learning_ , 8(3):279–292, 1992.
  * Alnæs et al. [2015] Martin S. Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E. Rognes, and Garth N. Wells.  The fenics project version 1.5.  _Archive of Numerical Software_ , 3(100), 2015.  doi: 10.11588/ans.2015.100.20553.
  * Hughes et al. [1986] Thomas JR Hughes, Leopoldo P Franca, and Marc Balestra.  A new finite element formulation for computational fluid dynamics: V. circumventing the babuška-brezzi condition: A stable petrov-galerkin formulation of the stokes problem accommodating equal-order interpolations.  _Computer Methods in Applied Mechanics and Engineering_ , 59(1):85–99, 1986.
  * Hughes and Franca [1987] Thomas JR Hughes and Leopoldo P Franca.  A new finite element formulation for computational fluid dynamics: Vii. the stokes problem with various well-posed boundary conditions: symmetric formulations that converge for all velocity/pressure spaces.  _Computer methods in applied mechanics and engineering_ , 65(1):85–96, 1987.
  * Codina [2002] Ramon Codina.  Stabilized finite element approximation of transient incompressible flows using orthogonal subscales.  _Computer methods in applied mechanics and engineering_ , 191(39-40):4295–4321, 2002.
  * Earl et al. [2007] David J Earl, CM Pooley, JF Ryder, Irene Bredberg, and JM Yeomans.  Modeling microscopic swimmers at low reynolds number.  _The Journal of Chemical Physics_ , 126(6):02B603, 2007.
  * Clift et al. [1978] Roland Clift, John R Grace, and Martin E Weber.  Bubbles, drops and particles.  _Academic Press_ , 1978.
  * Karp-Boss et al. [1996] L Karp-Boss, E Boss, PA Jumars, et al.  Nutrient fluxes to planktonic osmotrophs in the presence of fluid motion.  _Oceanography and Marine Biology_ , 34:71–108, 1996.
  * Fan et al. [2020] Jianqing Fan, Zhaoran Wang, Yuchen Xie, and Zhuoran Yang.  A theoretical analysis of deep q-learning.  In _Learning for Dynamics and Control_ , pages 486–489. PMLR, 2020\.