Stanislas Thiriet-Rupert

Jérôme Josse

David Perez-Pascual

Jason Tasse

Camille Andre

Lélia Abad

David Lebeaux

Jean-Marc Ghigo

Frédéric Laurent

Christophe Beloin
email: christophe.beloin@pasteur.fr

Analysis of in-patient evolution of Escherichia coli reveals potential links to relapse of bone and joint infections Running title (40 characters): in-patient evolution of E. coli in BJIs

Keywords: BJI, adaptation, relapse, Escherichia coli

40-word-or-less summary: Here, we identify that three E. coli strains involved in bone and joint infection (BJI) relapses adapted in-patient through various molecular mechanisms allowing to outcompete other bacteria or evading immune system. Understanding in-patient adaptation may help fighting relapses in BJIs.

Background

Bone and joint infections (BJIs), including prosthetic joint infections (PJI), are severe and difficult to treat infections affecting a growing number of patients, and showing 10 to 20% of relapse [1][2][3][4][5][6]. Although staphylococci are responsible for more than 50% of all BJI, 20% of them are caused by Enterobacterales, especially in the case of early and late acute infections [7], whereas Pseudomonas aeruginosa infections are associated with particularly difficult to treat BJI [6]. These data therefore support the need to study BJI caused by both Gram-positive and Gram-negative bacteria and for a better understanding of the in vivo adaptation leading to relapses.

During host colonization, pathogenic bacteria must face the endogenous microbiota or other co-infecting bacteria. This competition could involve toxin production [8,9] and other adaptive strategies to access resources [10] and evade the host immune system [11,12], eventually leading to a decrease virulence accounting for bacterial persistence in the host.

The genetic modifications underlying these adaptations are in general quite restricted upon natural host purifying selection [13,14]. In some extreme cases, mutator genotypes can evolve, thus facilitating within-host adaptation [15,16]. Extensive genome reductions and rearrangements [17][18][19] are also observed leading to the loss of non-essential genes in the host environment. However, there is an increasing body of evidence that these gene losses could also fuel adaptative evolution and microbial adaptation within patient is currently under intense scrutiny [20,21]. However, acquiring longitudinal samples and studying pathogens adaptation from the same patient in the context of infection relapses is difficult [START_REF] Trouillet-Assant | Adaptive processes of Staphylococcus aureus isolates during the progression from acute to chronic bone and joint infections in patients[END_REF]. In this study, we analysed three pairs of E. coli strains isolated from initial and recurrent BJIs in 3 patients to investigate how in-patient evolution occurred in these situations and provide insights on the potential mechanisms involved in infection relapse.

Methods

Selection of clinical strains

Initial/Relapse pairs of clinical E. coli strains used in this study were selected among the routine collections of the Bacteriology Department of Hôpital de la Croix-Rousse, Hospices Civils de Lyon and the Bacteriology laboratory at Cerballiance, Lyon. The first search criterion corresponded to collect E. coli strains isolated from BJI from 2017 to 2020. Then, initial/relapse pairs were formed when E. coli strain was isolated in two different samples at two different times in the same patient (initial and relapse) and at the same site of infection. Strain characteristics are available in Table 1. Antimicrobial susceptibility tests were performed by using the disc diffusion method (Supplementary Table 1).

Bacterial strains and growth conditions

Bacterial strains used in this study are listed in Supplementary Table 2. Bacteria were grown in Miller's Lysogeny Broth (LB) (Corning) or M63B1 minimum medium supplemented with appropriate antibiotic and carbon source as specified. Liquid cultures were incubated at 37°C or 30°C with 180rpm shaking. Solid cultures were done on LB with 1.5% agar supplemented with the appropriate antibiotics. Chemicals and media were purchased from Sigma-Aldrich.

Strain construction

Insertion of a Zeo-GFP and Zeo-mars genetic cassettes at the permissive lambda att attachment site as well as the deletion of fliC in I1 and R1 strains and the deletion of two adhesin coding genes in I2 strain were done using pKOBApra plasmid [START_REF] Chaveroche | A rapid method for efficient gene replacement in the filamentous fungus Aspergillus nidulans[END_REF] and lambda-red recombination.

Primers used for genetic constructions are listed in Supplementary Table 3.

CFU counting on biomaterials

The biofilm forming capacity was assessed on different orthopaedic biomaterials [stainless steel (AISI 316L), titanium (TA6V), medical grade silicone and ultra-high molecular weight polyethylene (UHMWPE)]. Uniform cylindrical pegs of orthopaedic biomaterial were manufactured by a supplier of medical and surgical orthopaedics equipment (Groupe Lépine, Genay, France) under the same conditions of manufacture of orthopaedic biomaterial and sterilised by gamma rays. Biomaterial pegs were fixed into the lid of 24-well flat bottom cell culture plates (Corning Inc., Corning, NY, USA). Bacterial biofilms were formed by immersing the lid-containing pegs into 2mL of cultures at OD600=0.05 and incubated in LB for 24 hours at 37°C in humid atmosphere. After 3 rinsing steps in PBS, the pegs were dropped in a Falcon tube containing 5mL of PBS for biofilm disruption using sonication for 5 minutes (BACTOSONIC 14.2, BANDELIN, Berlin, Germany). Then, viable cell count was determined by serial dilution and plating.

Bacterial adhesion to MG63 human osteoblastic cells

The osteoblastic cell line MG63 (LGC standards, Molsheim, France) was used to test bacterial adhesion. Cells were cultured in 75cm² flasks (T75, BD Falcon, Le Pont de Claix, France) at 37°C under 5% CO2, in a culture medium composed of DMEM (Dulbecco's Modified Eagle Medium containing D-glucose, L-glutamine, pyruvate) supplemented with 10% fetal calf serum (FCS), penicillin (100µg/mL) and streptomycin (100µg/mL), all from Gibco (Paisley, UK).

Cells were passaged once a week and used until passage 25 at the most. The day before the infection, MG63 cells were seeded at 100,000 cells per well in a 24-well plate and the different bacterial strains were inoculated in LB. The day of the infection, overnight bacterial cultures were centrifugated, resuspended and adjusted in DMEM + 10% FCS to reach a MOI (multiplicity of infection) of 100:1 and added to cell culture. After 2 hours of infection, cells were washed twice with PBS and incubated and lysed by osmotic shock in sterile water and lysates were plated on tryptone soy agar plates (Biomérieux) using the Easy Spiral® automaton (Interscience, Saint-Nom-la-Bretèche, France). CFU were counted after 18 hours of incubation at 37°C and results were presented as CFU for 100 000 cells.

Survival to human serum

For each strain, 50µL at OD600=0.4 was mixed in a well of a 96-well plate with 50µL of active human serum. The same process was followed for a control with heat inactivated serum (HIS) (30 minutes at 55°C). Plate was then incubated at 37°C under shaking (450rpm) for 1h30. Each well was then serially diluted and plated on LB for CFU count and the survival of each strain was expressed as the ratio between the CFUs in active serum divided by the CFUs in HIS.

Zebrafish infection

Overnight bacterial cultures were centrifuged, washed, and resuspended at the desired concentration in PBS. Anesthetized zebrafish larvae were microinjected intravenously 52-58 hours post-fertilization with 1nL of each bacterial suspension at a final concentration of 2x10 3 CFU [START_REF] Colucci-Guyon | Strategies of professional phagocytes in vivo: unlike macrophages, neutrophils engulf only surface-associated microbes[END_REF]. Infected larvae were transferred into individual wells containing 1mL of Volvic water (Danone, France) and 0.003% 1-phenyl 2-thiourea (used to improve optical transparency in zebrafish) in 24-well plates, incubated at 28°C, and regularly observed under a stereomicroscope. Protocols describing zebrafish maintenance and analysis, see the corresponding supplementary method section.

Competition assay for biofilm formation

The pair of strains used in competition assay carried two different fluorescent tags and were inoculated in LB medium supplemented with appropriate antibiotic and incubated overnight at 37°C under shaking. Cultures were then adjusted to OD600 = 0.05 and mixed in 1:1 ratio (verified by FACS). Each competition mix was then inoculated in a PVC plate and incubated at 37°C for 16 hours. The supernatant was removed from each well, the biofilm resuspended in 100µL of PBS and the proportion of each strain was assessed by FACS and the relative fitness calculated as follows:

𝑟𝑒𝑙𝑎𝑡𝑖𝑣𝑒 𝑓𝑖𝑡𝑛𝑒𝑠𝑠 = - 𝑇𝑎𝑟𝑔𝑒𝑡 !"#$ 𝑇𝑎𝑟𝑔𝑒𝑡 !%& 0 - 𝑟𝑒𝑓 !"#$ 𝑟𝑒𝑓 !%& 0
, where Targetpre and Targetpost are the cell concentrations of the target strain in the mixed culture before and after the overnight incubation, and refpre and refpost are the cell concentrations of the corresponding reference strain. In order to account for growth differences due to the expression of one or the other fluorescent tag (supplementary Figure 1), competitions were run with both possible combinations of fluorescent tags. This protocol was used to compete R1 and R3 strains against their respective ancestors. However, because GFP and RFP tags were not expressed in I2 and R2 strains, we used an alternative strategy based on PCR amplification and Sanger sequencing of a SNP locus allowing to destinguich both strains (see supplementary method section).

Competition assay for growth in liquid cultures

The strains were incubated at 37°C in the appropriate medium and adjusted to OD600=0.05. A 1:1 ratio mixed culture was verified by FACS and incubated overnight at 37°C with shaking.

The concentrations of target and reference strains was measured using FACS and relative fitness were calculated as above.

SNP analysis

Mutations in the adapted strain were identified by using the genome assembly of its ancestral strain (see supplementary methods) as a reference. breseq version 0.35.0 [START_REF] Deatherage | Identification of mutations in laboratory evolved microbes from next-generation sequencing data using breseq[END_REF] with the consensus mode was used with default parameters. As a confirmation of this analysis, Freebayes version 1.3.3 [START_REF] Garrison | Haplotype-based variant detection from short-read sequencing[END_REF] was also used to identify SNPs using the following stringent thresholds: minimum depth of 20 reads, a maximum depth corresponding to the average depth over the genome plus one standard deviation, at lease 10 forward and 10 reverse reads supporting the SNP, minimum variant quality of 30, minimum base quality of 20. Only variants identified in both analyses were validated. The resulting identified mutations were filtered out if they were too close to each other (less than 50bp). These mutations are usually the result of misaligned reads, often due to repetitive regions [START_REF] Ghalayini | Evolution of a Dominant Natural Isolate of Escherichia coli in the Human Gut over the Course of a Year Suggests a Neutral Evolution with Reduced Effective Population Size[END_REF]. Because these regions are not 100% identical, they are disturbing the mapping process resulting in false positives. Mutations found in mobile elements were kept if only one copy was present in the genome.

Statistical analysis

Figures and statistical tests were performed using either Prism 9.5.1 for Mac OS X (GraphPad Software, Inc.) or Rstudio with R version 4.2.1 using the R package GGplot2 [START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF].

Methods related to the supplementary results are available in supplementary methods.

RESULTS

I-In-patient evolved Escherichia coli clones showed no antibiotic resistance acquisition but signs of natural selection.

In order to investigate potential cases of in-patient adaptation, we obtained three pairs of E. coli strains sampled in three patients presenting relapsing BJIs. Each pair is composed of a sample at the time of the initial infection (strains hereafter denominated as I1 to I3 for cases 1 to 3) and a second sample corresponding to the relapse of this infection (strains hereafter denominated as R1 to R3 for cases 1 to 3). The genome of each strain was sequenced and analyzed to confirm that the infection and relapse strains were indeed the same and not a case of re-infection (Table 1, see supplementary methods). Interestingly, antimicrobial susceptibility testing (supplementary Table 1) and in silico analysis (supplementary Data 1) did not reveal any acquisition of antibiotic resistance associated with the relapse events. Consequently, each strain pair was examined for other phenotypic and genetic features that could have been selected

for to allow the infection relapses.

The first notable difference identified by genomic analysis was the number of gained and lost genes (supplementary Table 4), some of which corresponding to potential plasmids (supplementary Data 2, see supplementary methods). Particularly, the R2 strain that displayed the longest in-patient evolution underwent an important genome reduction with the loss of 3.4% of its size corresponding to 170 coding sequences involved in various functions such as transport, iron acquisition, carbohydrate acquisition and utilization, defense against competitors with Microcin H47 production-related genes and three adhesins (supplementary Data 2).

The number of mutation per base per year was quite high for all three relapse strains (1.11x10 - 10 the gut, [START_REF] Ghalayini | Evolution of a Dominant Natural Isolate of Escherichia coli in the Human Gut over the Course of a Year Suggests a Neutral Evolution with Reduced Effective Population Size[END_REF]). In addition, even if no functional convergence was observed in the mutated genes, the predominance of loss of function and non-synonymous mutations suggests that positive selection was at play in these cases (Table 2, see supplementary methods).

II-In vivo adaptation did not modify biofilm formation capacities.

Because strain motility could impact the virulence and biofilm forming capacity, we tested each strain on low-agar LB plates at 30°C and 37°C, which showed that the second and third pairs of strains were non-motile (supplemental Figure 2, see supplementary methods).

Both I1 and R1 strains display temperature dependent motility very likely due to a flagellin production restricted to 30°C (supplementary Figure 2 and3, see Supplementary Results).

Biofilm formation was assessed on various surfaces showing that the R1 relapse strain formed 1.71 more biofilms than I1 strain when grown on silicon coupons (Figure 1.A). No difference was identified for the I2-R2 pair (Figure 1.B), while the I3 strain showed a 1.49 increase in biofilm formation on silicon (Figure 1.C). However, in all cases the tested clinical strains showed a weak biofilm forming capacity, all of them being equivalently or even less efficient than E. coli K12 strain MG1655 used as control for weak biofilm formation.

Interestingly, R1 strain showed a slight decrease of its adhesion capacity to osteoblasts (Figure 1.D). However, no evident genetic modification could be linked to this phenotype. No difference was observed either in term of auto-aggregation or production of curli amyloids and cellulose (supplementary Figure 4, see supplementary methods and Supplementary Results).

Taken together, these results show that very limited differences in biofilm formation were detected within each pair suggesting that a modification of biofilm formation and biofilmrelated phenotypes was unlikely to be involved in in vivo adaptation in these three studied cases.

III-Loss of adhesins induces serum resistance in R2 strain.

To assess the strain capacity to challenge to the host immune defenses, we measured the survival of all strains to exposure to active human serum. While we observed no difference between I3 and R3 strains, the relapse strain R1 had a slight decrease in serum resistance compared to I1 strain and the relapse strain R2 had a clear increase in serum resistance relatively to its ancestral strain I2 (Figure 2.A). Although we identified no difference in the I2 vs R2 lipopolysaccharide profile (outer-membrane component targeted by the complement) (Figure 2.B), the loss of two adhesins in R2 could be linked to this phenotype (supplementary Data 2) [11]. Indeed, deletion of each adhesin in I2 strain led to reduced serum resistance suggesting that the loss of both adhesins could have been selected through the conferred immune system evasion, contributing to the strain adaptation to the host environment (Figure 2.C).

Unfortunately, we were unable to assess the additive effect of these mutations since we could not reconstruct the double mutant in the I2 strain despite repeated attempts.

IV-Strain virulence was not targeted during in-patient selection.

To gain insight into the potential evolution of the virulence of the studied strains, we assessed their cytotoxicity on HeLa cells as well as hemolysis capacities and showed no differences, except for a slightly lower hemolysis induced by R2 strain after 24h of incubation in presence of defibrinated horse blood (Figure 3.A.B).

Potential differences in in vivo virulence were tested in a zebrafish (Danio rerio) larval model [START_REF] Wiles | Use of Zebrafish to Probe the Divergent Virulence Potentials and Toxin Requirements of Extraintestinal Pathogenic Escherichia coli[END_REF]. While no differences were found in bacterial load for all pairs (Sup Figure S5.A.B.C), zebrafish larvae exposed to R2 relapse strain had a slightly lower survival, although nonsignificant, 24 hours post infection (hpi) than with its related ancestor I2, even though both strains ultimately induced the same level of mortality 48 hpi (Figure 3.C.D). No difference was identified concerning the zebrafish innate immune response level (Figure 3.E). These results suggest that the R2 relapse strain induced a slightly quicker mortality, which does not seem to be linked to a difference in colonization ability.

V-Growth optimization and colicin production as strategies to outcompete other bacteria.

To assess if a difference in growth capacity could advantage the strains isolated after relapses, they were pitted against their related ancestor in 1:1 ratio competition experiments in biofilm and liquid conditions. First, no difference in fitness was observed in biofilm condition, which further straightened the hypothesis that biofilm formation was not targeted during the in vivo adaptation. In liquid, while R1 strain slightly outcompeted its ancestor I1, R2 strain had a striking fitness advantage over I2 (Figure 4.A).

To further investigate this phenotype, we competed I2 and R2 strains in presence of different sugars as well as of an iron chelator (2,2'-Bipyridyl) based on identified genetic modifications and especially gene losses related to uptake and metabolism of these substrates (supplementary Data 2). Surprisingly, despite these gene losses, R2 strain showed no growth defect and even globally had an advantage through a capacity to better sustain growth during stationary phase (supplementary Figure 6), as confirmed by competition assays (Figure 4.B).

Therefore, extensive gene loss seems to have resulted in a global growth optimization, which could have helped the strain to outcompete other bacteria including its ancestor.

Another strategy to outcompete its kin and promote subsequent adaptation was identified in the relapse strain R1, which acquired a colicin-producing plasmid allowing the production of colicin E1. An overlay assay indeed confirmed the acquired capacity of R1 relapse strain to inhibit I1 strain growth (supplementary Figure 7).

Discussion

Treatment failure is a common issue of BJIs, allowing bacteria to remain and adapt to the host, ultimately increasing the risks of relapse and the difficulty to cure the patient. A better understanding of the mechanisms underlying pathogens adaptation to the host environment would allow improving the management of relapsing BJIs. Here we investigated three cases of BJIs relapse involving pathogenic E. coli.

Our genetic and phenotypic analysis revealed diverse functions arising from bacterial adaptation to the host environment during BJIs. We found that adapted strains isolated from relapses displayed a relatively high mutation rate, which is a known adaptive advantage increasing the advantageous mutation probability that could be subsequently selected [13,[START_REF] Giraud | Costs and benefits of high mutation rates: adaptive evolution of bacteria in the mouse gut[END_REF][START_REF] Tenaillon | Tempo and mode of genome evolution in a 50,000-generation experiment[END_REF]. However, no mutations were found in genes associated with DNA damage repair, suggesting that these high mutation rates were intrinsic to these strains and could have further promoted the selection of beneficial mutations, partly explaining their adaptive success.

Although biofilms are known to be frequently involved in infections and relapses, none of the three investigated cases showed an increased biofilm formation capacity that could have played a role in the relapses and all strains were even characterized as poor in vitro biofilm formers.

While a strong in vitro biofilm formation is often linked to the severity of clinical outcomes in prosthesis-related infections [START_REF] Malchau | Biofilm properties in relation to treatment outcome in patients with first-time periprosthetic hip or knee joint infection[END_REF][START_REF] Sanchez | Biofilm formation by clinical isolates and the implications in chronic infections[END_REF][START_REF] Post | Comparative Genomics Study of Staphylococcus epidermidis Isolates from Orthopedic-Device-Related Infections Correlated with Patient Outcome[END_REF], it did not seem to be a critical factor in this study.

The host immune system exerts a critical selection pressure on pathogenic bacteria. In I1 and R1 strains, the inhibition of flagellin production at 37°C could minimize their recognition by the host's immune defenses [START_REF] Smith | Toll-like receptor 5 recognizes a conserved site on flagellin required for protofilament formation and bacterial motility[END_REF][START_REF] Bambou | In vitro and ex vivo activation of the TLR5 signaling pathway in intestinal epithelial cells by a commensal Escherichia coli strain[END_REF]. This could not only contribute to immune system evasion but also to save energy through a down-regulation or suppression of flagellin production, as

shown to be the case in the context of gut infection by E. coli [START_REF] Giraud | Dissecting the genetic components of adaptation of Escherichia coli to the mouse gut[END_REF][START_REF] Leatham | Mouse intestine selects nonmotile flhDC mutants of Escherichia coli MG1655 with increased colonizing ability and better utilization of carbon sources[END_REF].

This immune system evasion strategy is even clearer in the case of the R2 relapse strain, which showed increased survival to human serum through the loss of two adhesins. This suggests that in these cases, evading the host's immune system may be more critical than maintaining virulence as shown to happen during in vivo adaptation [START_REF] Winstanley | Pseudomonas aeruginosa Evolutionary Adaptation and Diversification in Cystic Fibrosis Chronic Lung Infections[END_REF][START_REF] Viberg | Within-Host Evolution of Burkholderia pseudomallei during Chronic Infection of Seven Australasian Cystic Fibrosis Patients[END_REF]. This could allow pathogens to stay longer in the host as shown in the infection experiments performed in zebrafish, where no difference was observed.

In addition, the R2 strain presented an interesting case of gene loss. Genome reduction is usually considered as a strategy selected because energy consumption is reduced via the loss of genes that are not used or under selective pressure [START_REF] Moran | Microbial minimalism: genome reduction in bacterial pathogens[END_REF][START_REF] Weinert | Why Might Bacterial Pathogens Have Small Genomes?[END_REF]. In the case of the R2 strain, numerous gene losses seem to result in a global growth capacities optimization, potentially associated to lower the cost for gene production and lower amount of DNA to replicate. It is possible that some nutritional resources could be directly acquired from the new host environment instead of being synthesized by the cells and that the loss of genes involved in this production optimizes the cell metabolism [START_REF] Souza | Experimental Evolution of Metabolic Dependency in Bacteria[END_REF]. Gain of growth capacity in the context of an infection could be of paramount importance to outcompete other endogenous or pathogenic bacteria and thrive in the host environment.

The polymicrobial context of some of these BJIs is also an important parameter to stress since the interactions between the different bacterial species or strains could play an important role in the relapses and the mechanisms shaping bacterial adaptation to the host environment [START_REF] Tay | Polymicrobial-Host Interactions during Infection[END_REF][START_REF] Boldock | Human skin commensals augment Staphylococcus aureus pathogenesis[END_REF]. This study may therefore have missed some adaptation that rely on bacteria-bacteria direct or indirect interactions. In addition, as for any longitudinal study, there is always a possibility that the strains causing the very initial infection already evolved at the time of sampling when infection was clinically detected, and some of these mutations, that would be undetected in our study, could potentially contribute to the success of the relapse strains.

Finally, the small number of strains studied, and the diversity of BJI cases, make it challenging to draw definitive conclusions about typical strain profiles and mechanisms involved in in vivo adaptation. However, it seems that a higher mutation rate allowing a better adaptability as well as evolving immune system evasion and outcompeting other bacteria are relevant strategies selected in patient. Further investigations with a larger cohort would help determine whether these diverse adaptive strategies constitute a more general rule of E. coli pathogenic strain in BJIs.

ETICS STATEMENT

Animal experiments were performed according to European Union guidelines for handling of laboratory animals (http://ec.europa.eu/environment/chemicals/lab_animals/home_en.htm) and were approved by the Institut Pasteur Animal Care and Use Committee and the French Ministry of Research (APAFIS#31827). to its ancestor (n=3 biological replicates). Values higher than 1 show higher hemolysis capacity for the relapse strain, while values lower than 1 show lower hemolysis capacity for the relapse strain. Statistical significance according to non-parametric two-tailed Mann-Whitney test. * p<0.05. (C) Monitoring of zebrafish survival (no statistical significance was detected according to log-rank test), (D) bacterial burden into the zebrafish (no statistical significance was detected according to non-parametric two-tailed Mann-Whitney test) and (E) in vivo fluorescence of macrophages and neutrophils for zebrafish infected by I2 and R2 strains as well as non-injected supernatant from a culture in presence of Mitomycin C to induce colicin production. In both cases, LB is used as a negative control. Supplementary Table 1. Antibiotic susceptibility profiles. Resistant (R), sensitive (S) or intermediate (I) status of each strain to commonly used antibiotics assessed by disc diffusion assay and following EUCAST breakpoints recommendations.

Antibiotic I1 R1 I2 R2 I3 R3 Ampicillin R R R R R R Cefoxitin S S S S R R Meropenem S S S S S S Ticarcillin R R R R R R Ticarcillin + clavulanic acid S S R R R R Cefepime S S R R R R Ertapenem S S S S S S Mecillinam S S S S S S Cefotaxime S S R R R R Amoxicillin + clavulanic acid S S R R R R Ceftazidime S S I I R R Piperacillin R R R R R R Tazocillin (Piperacillin + Tazobactam) S S R R R R Aztreonam S S R R R R Imipenem S S S S S S Cefalexin R R S S S S Tobramycin S S S S R R Amikacin S S S S S S Norfloxacin R R S S R R Moxifloxacin R R S S R R Ciprofloxacin R R S S R R Co-trimoxazole (Trimethoprim + Sulfamethoxazole) R R S S R

Supplementary Methods

RDAR morphotypes on congo red plates

Two microliters of an overnight culture grown at 37°C in LB medium under shaking (with appropriate antibiotics when needed) was spotted onto LB plates (without NaCl) containing 0.004% Congo Red and 0.002% brilliant blue. The spotted drops were allowed to dry, and the plates were incubated for 24 and 48 hours at 30 or 37°C.

Auto-aggregation curves

Cultures were grown as described above, and 1mL aliquots were diluted to an OD600=3 in spent medium in order to prevent growth during experiment. Once diluted, all the tubes were vortexed and an aliquot of 50µL was removed, 1cm below the top of the culture, and diluted with 50µL of LB prior to OD600 measurement. The bacteria were then left to aggregate at RT for 6 hours, while repeating the measurement of OD600 1 cm below the top of the culture every hour.

Motility assay

Bacterial strains were incubated overnight in LB at the specified temperature with shaking.

Then, 2µL of overnight culture was inoculated in swim plates (1% tryptone, 0.25% NaCl, 0.25% Eiken agar) for 24 hours at 37°C or 30°C. After incubation, motility was assessed by checking the spread of inoculated bacteria into the motility agar.

Biofilm microtiter plate CV assay

Bacteria were grown overnight and diluted to an OD600=0.05 and 100µL were inoculated in technical triplicates in either a polystyrene Greiner 96-well plate or a polyvinyl chloride (PVC) round bottom 96-well plates (Corning). After 24 hours at 37°C, the biofilms were washed twice by slowly pipetting PBS in each well to remove non-attached bacteria. Biofilms were then resuspended in 100µL of PBS by vigorous pipetting and viable cell count was determined by serial dilution and plating.

Growth capacity

An overnight culture of each strain was diluted to OD600 of 0.05 in appropriate medium. One hundred and fifty μL aliquots were inoculated in a 96-well plate. The plates were then incubated in a TECAN Infinite M200 Pro spectrophotometer (Männedorf, Switzerland) for 18 hours at 37°C with shaking of 2mm amplitude. The absorbance of each culture at 600 nm was measured every 15 minutes. Growth curve of each strain was measured using 3 biological replicates.

The R package GrowthCurver [1] was used to infer growth-related parameters for all strains based on the growth curves.

Cytotoxicity on HeLa cells

The cytotoxic activity on Hela cells of each bacterial strain was analyzed using the CellTiter-Blue® Cell Viability Assay (Promega) following manufacturer instructions. HeLa cells were cultured in complete DMEM (Life Technologies) supplemented with 10% FCS, 100µg/mL penicillin/streptomycin, and 10mM L-glutamine at 37C in the presence of 5% CO2. For infection assay, the HeLa cells were incubated in 96-well plates to reach a density of 10 5 cells.

Each bacterial strain was grown overnight in LB as described above and then washed twice in PBS, to finally adjust the bacterial density to 10 7 CFU/ml. These bacterial suspensions were used to infect HeLa cells by adding 100 µl of each to the wells (MOI of 1:100). The plates were incubated at 37C in the presence of 5% CO2 for 2 hours. Non-infected cells and the E. coli K12 MG1655 strain were used as negative controls. At the end of the incubation time, the wells were washed twice with PBS to remove the bacteria. Then, 100 µl of PBS and 20 µl of CellTiter Blue solution were added to each well and incubated for 2 hours at 37C and 5% CO2. The cell viability was measured by recording the fluorescence at 560/590 nm using a Tecan Infinite-M200-Pro spectrophotometer.

Hemolysis capacity

The hemolytic activity of each strain was assessed by mixing defibrinated horse blood (bioMerieux, Marcy l'Etoile, France) diluted to 10% (v/v) in PBS with a total of 10 8 bacteria in PBS. After 2, 6 and 24 hours of incubation at 37°C with shacking, the non-lysed erythrocytes were pelleted by centrifugation at 1500g for 2 minutes, and the OD450 of 150μL of the cell-free supernatants was measured. The hemolytic capacity of each adapted strain was expressed relatively to its ancestral strain.

Luria-Delbrück fluctuation test

The mutation rate of I1, R1 and MG1655 strains were assessed using Luria-Delbrück fluctuation test [2]. For each strain, 60 wells of a 96-well plate were inoculated with ~100 bacteria in a volume of 200µL in LB before incubating the plate in a humid atmosphere for 16h at 37°C under 450 rpm shaking. After the cultures had reach saturation, ten randomly chosen wells were sampled to evaluate the number of total bacteria per well by serial dilution and plating. The total volume of each of the remaining 50 wells was plated on individual LB plates supplemented with 2xMICs of Rifampicin to evaluate the number of resistant cells. The colony count was then used to assess the mutation rate by a maximum likelihood approach available on FALCOR web tool (https://lianglab.brocku.ca/FALCOR/) [3].

Zebrafish care and maintenance

Homozygous Tg(mfap4::mCherryF) (ump6Tg) [4] Tg(mpx::GFP) i114 [5] double transgenic fishes were raised in Institut Pasteur zebrafish facility. Eggs were obtained by natural spawning, bleached according to standard protocols, and then kept in Petri dishes containing Volvic spring water and, from 24 hpf onwards, 0.003% 1-phenyl-2-thiourea (PTU) (Sigma-Aldrich) was added to prevent pigmentation. Embryos as well as infected larvae were reared at 28°C. Larvae were anesthetized with 200 μg/mL of buffered tricaine methane sulfonate (MS-222, Sigma-Aldrich) during the injection procedure as well as during in vivo imaging and processing for bacterial burden evaluation.

Evaluation of the bacterial burden in infected zebrafish larvae

Infected zebrafish larvae were collected at 0, 24, and 48 hpi and lysed for the evaluation of the bacterial burden as previously described [6]. Each larva was placed in an individual 1.5 mL Eppendorf tube and anesthetized with tricaine (200 μg/mL), washed with 1 mL of sterile water, and placed in 150 µL of sterile water. Larvae were then homogenized using a pestle motor mixer (Argos). Each sample was transferred to an individual well of a 96-well plate and 10 X serial dilutions were performed. For CFU enumeration and to assess plasmid stability throughout the infection kinetics, serial dilutions of lysates were plated on MacConkey agar plates supplemented with lactose plates and incubated overnight at 37°C.

Visualization of macrophages and neutrophils in infected zebrafish larvae

Visualization of macrophages and neutrophils on living transgenic reporter larvae was performed upon infection as we previously described [7]. Briefly, bright field, GFP, and RFP images of whole living anesthetized larvae were taken using a Leica Macrofluo Z16 APOA (zoom 16:1) macroscope equipped with a Leica PlanApo 2.0 X lens, and a Photometrics CoolSNAP HQ2 camera. Images were captured using the Metavue software version 7.5.6.0 (MDS Analytical Technologies). After capture of images, larvae were washed and transferred to a new 24-well plate filled with 1mL of fresh water per well and incubated at 28 °C.

Overlay assay colicin

I1 and R1 strains were incubated overnight at 37°C with shaking in LB supplemented with 40 µg/mL of Mitomycin C to induce colicin production. After incubation an overlay of each strain was realized by pouring 6mL of overnight culture diluted at OD600=0.1 onto a LB agar plate and removing it after 5 minutes. Plates were allowed to dry and a 10µL drop of supernatant from each strain culture filtered using 0.22uM filter was spotted on each overlay. After an additional drying step, the plates were incubated overnight at 37°C and the production of colicin was identified by the presence of an inhibition halo when the filtered supernatant was spotted.

LB was spotted as a negative control and filtered supernatant from both strains was also spotted on overlay of each of them to be sur that the growth inhibition halo was not due to the presence of Mitomycin C in the filtered supernatant.

Whole genome sequencing and analysis

Prior to genome extraction, a single colony of each strain was inoculated in LB medium for over-day culture till the OD600 reached around 1.0 (ca. 5.0x10 8 bacteria/mL). The bacterial cells were collected from 2mL of the culture and the genomic DNA was extracted using the Wizard Genomic DNA Purification Kit (Promega). Samples were then sent to Novogene Co., Ltd to be sequenced using Illumina NovaSeq 6000. Sequencing reads were pre-processed to remove lowquality or artefactual bases. We used fqCleaner v.0.5.0, a mini workflow implemented in Galaxy [8] (https://galaxy.pasteur.fr/) to process fastq files (quality trimming, duplicate and artifact filters). Draft genomes were de novo assembled using SPAdes version 3.15.0 [9] and annotated using Proka version 1.14.5 [10]. Assembly quality and completeness was assessed using BUSCO version v5.4.2 [11]. The assembled genomes were used for in silico typing and characterization of each strain using tools available at the Centre of Genomic Epidemiology (http://www.genomicepidemiology.org/): Resfinder [12][13][14], MLST [14,15], SerotypeFinder [16], CHtyper [17]. Phylogroups were identified in silico using ClermonTyping tool [18]. In order to identify contigs corresponding to potential plasmids, each contig was used as a query in a BLASTn search using blast+ version 2.2.31 [14] against the COMPASS database [19].

Only hits with 95% identity and 80% coverage of both the query and the hit sequence were kept.

Alternative strategy for competition assays of I2 vs R2 strains

Even if we managed to introduce GFP and RFP (mars) tags in I2 and R2 strains, the related proteins were likely not expressed since thefluorescence was impossible to detect.

Consequently we used an alternative strategy to assess the proportion of each strain. We took andvantage of the fact that both strains can be differentiated by SNPs. An aliquot of the mixed culture containing 10 7 cells was pelleted and resuspended in 50µL of sterile distiled water and boiled for 10 minutes. Specific primers were then used to amplify a SNP locus in butC gene and Sanger sequence the resulting DNA fragment. Then, using QSVanalyzer [20], we calculated the frequeny of each mutation and used it to compute the relative fitness as above.

To assess the accuracy of this strategy, we used strains I1 and R1 taged with different fluorescent marker to compare strains proportions as assessed by FACS and QSVanalyzer and found a very good correspondance between both methods (Supplemenatry Figure S6), therefore validating our alternative approach.

Western blot analysis and LPS gel

The heat extracted proteins from 1mL at OD600=1 culture were suspended in 1×Laemmli buffer and incubated for 5 min at 95°C. The protein extracts were run on Mini-PROTEAN TGX Stain-Free precast Gels (BioRad) in 1×TGX buffer and then transferred to nitrocellulose membrane using a Trans-Blot ® Turbo Transfer System (BioRad). Blocking was performed in a 5% solution of dry milk and 0.05% Tween 1×PBS (1×PBST) overnight at 4°C with agitation. The membranes were then incubated in 1×PBST with either rabbit anti-flagellin (kindly given by Dr Monica Rolando) or mouse anti-RNAPa (BioLegend #663104) diluted 1/5,000 for 1h at room temperature with agitation. Membranes were washed in 1×PBST and then incubated with the secondary antibody (anti-rabbit or anti-mouse IgG conjugated with horse radish peroxidase at 1:10 000, Promega). The membranes were then revealed using an ECL kit (GE Healthcare) and the iBright TM CL1500 system (Thermofisher).

LPS analysis on Mini-PROTEAN TGX Stain-Free precast Gels (BioRad) were realized as in [21].

Figure 1 .

 1 Figure 1. Biofilm formation capacity for infection and relapse strains from case 1 (A), 2 (B) and 3 (C). Six different surfaces were tested, E. coli K12 MG1655 was used as a control for weak biofilm formation and E. coli TG1 as a control for strong biofilm formation (n=3 to 4 biological replicates). (D) Adhesion capacity on osteoblast human cells (n=3 biological replicates). For each panel, statistical significance was only assessed between the infection strain and its related relapse strain according to non-parametric two-tailed Mann-Whitney test. * p<0.05; ** p<0.01.

Figure 2 .

 2 Figure 2. Serum survival capacity. A. Normalized serum survival capacity of each strain expressed as percentage of the control for growth in Heat Inactivated Serum (HIS). E. coli strains CFT073 and MG1655 were used as positive and negative controls, respectively (n=3 biological replicates). Statistical significance was only assessed between the infection strain and its related relapse strain according to non-parametric two-tailed Mann-Whitney test. * p<0.05. B. LPS profile of I2 and R2 strains. E. coli MG1655 was used as a negative control and the same strain in which the O-Antigen was reconstructed was used as a positive control. C. Relative Serum survival of I2 ancestral strain, R2 relapse strain as well as two deletion mutants constructed in the I2 ancestral strain: one for the afa fimbrial operon and a second for the yadA-like encoding gene. Results are expressed relatively to R2 relapse strain (n=4 biological replicates). Statistical significance according to non-parametric two-tailed Mann-Whitney test. * p<0.05; ** p<0.01.

Figure 3 .

 3 Figure 3. Virulence of the strains. (A) Cytotoxicity test on HeLa cells. All six adapted and ancestral strains were tested as well as MG1655 strain as a negative control. The values for each strain were expressed relatively to the control without bacteria (n=3 biological replicates). No statistical significance was detected according to nonparametric two-tailed Mann-Whitney test. (B) Hemolysis capacity in liquid culture of each relapse strain relatively

  Figure 1

Table 1 .

 1 Cases and strain information. Table summarizing the type and location of infection as well as type of sample. Each strain was typed in silico and, when the infection was polymicrobial, the other species involved are listed. ANI: Average Nucleotide Identity, MLST: MultiLocus Sequencing Typing.

	Case N°	Strain	Infection status	Infection site	Sample type	Diagnosis	Time before relapse	Antibiotic susceptibility profile	Phylogroup	Serotype	ANI	MLST #1	MLST #2	fimH	fumC	Infection type	Other pathogens if polymicrobial
					Deep
		I1	Initial		tissue by			Penicillinase								Monomicrobial
					swab	Septic
	1			Hand	Deep	arthritis 4th finger (no	51 days		B1	O76:H28	99.99%	224	730	FimH3 9	FumC4		Enterococcus faecium,
		R1	relapse		tissue by	material)		Identical								Polymicrobial	Pseudomonas
					swab												aeruginosa,
																	Citrobacter freundii
	2	I2 R2	Initial relapse	Spine	Washing liquid biopsy	Vertebral osteomyelitis (no material)	748 days	ESBL-producing Identical	G	O24:H4	99.99%	117	48	FimH9 7	FumC4 5	Monomicrobial Polymicrobial	Cutibacterium acnes
	3	I3 R3	Initial relapse	Right hip	Bone Deep tissue	Prosthetic joint infection	66 days	ESBL-producing Identical	A	O89/O16 2:H9	99.99%	167	2	/	FumC1 1	Monomicrobial Monomicrobial

Table 2 .

 2 List of all strains used in this study. . Mossoro C, Glaziou P, Yassibanda S, et al. Chronic Diarrhea, Hemorrhagic Colitis, and Hemolytic-Uremic Syndrome Associated with HEp-2 Adherent Escherichia coli in Adults Infected with Human Immunodeficiency Virus in Bangui, Central African Republic. J Clin Microbiol. 2002; 40(8):3086-3088. 5. Da Re S, Ghigo J-M. A CsgD-independent pathway for cellulose production and biofilm formation in Escherichia coli. J Bacteriol. 2006; 188(8):3073-3087.

	R

1. Mobley HL, Green DM, Trifillis AL, et al. Pyelonephritogenic Escherichia coli and killing of cultured human renal proximal tubular epithelial cells: role of hemolysin in some strains. Infect Immun. 1990; 58(5):1281-1289. 2. Rendueles O, Beloin C, Latour-Lambert P, Ghigo J-M. A new biofilm-associated colicin with increased efficiency against biofilm bacteria. ISME J. 2014; 8(6):1275-1288. 3. Francius G, Polyakov P, Merlin J, et al. Bacterial Surface Appendages Strongly Impact Nanomechanical and Electrokinetic Properties of Escherichia coli Cells Subjected to Osmotic Stress. PLOS ONE. Public Library of Science; 2011; 6(5):e20066. 4

Table 3 .

 3 List of all primers used in this study.

	I1/R1 and I3/R3 lATT Zeo-mars and Zeo-GFP construction and verification
	lambda-ATT. B1.500-3	CCCTGATACTCACCAGGCATCAC
	Rev_K12-ZeoXFP-lATT	GTCACGCCAAAAGCCAATGC
	Fw_K12-ZeoXFP-lATT	GCAAGCGCCTCGATTACTGC
	lambda-ATT.A1.500-5	CGATGGCGATAATATTTCACC
	Sanger sequencing-based fitness assessment in I2/R2
	Fw_butC	CAGTCGGTTATTGCTGGCT
	Rv_butC	CACACTTTGAGGCGACATTG
	Sanger sequencing-based fitness assessment calibration vs FACS in I1/R1
	Fw_yehT	CAACATTTGATGCTGGAGATCG
	Rv_yehT	GCTTCATCAATTGGCTTCAGC
	afaBCDE deletion in I2 and verification
	Fw_afaDel_size_ZeoFRT	GATTCAGCACCTCAGTCAGAC
	Fw_Del_afaB-E_ZeoFRT	GGTTGCTATTAACATATTTATAAAAGGTTTATTTGCCTTCAGGATAAAATATGTAGGCTGGAGCTGCTTC
	Rv_Del_afaB-E_ZeoFRT	CAAAGGGAGCATATAGCCCCCTTCTTTCATCAGGTTAATTTTCCAGAATATCCTCCTTAGTTCCTATTCC
	Rv_afaDel_size_ZeoFRT	GCGAGGAAGATTTCTCTTGTAGG
	yadA-like deletion in I2 and verification
	Fw_Size_yadA_Del	TGCTGCTTAAATTGAGGTCGC
	Fw_Del_yadA_ZeoFRT	CTCCGCACTACCGTTCTGGCTGGTGAAGTAATAATACAGGAGAATAACAGATGTAGGCTGGAGCTGCTTC
	Rv_Del_yadA_ZeoFRT	TCCCGGGGGACTCCCCCCGGGACAGATATCTAAATCCTGACCTGGAATATCCTCCTTAGTTCCTATTCCG
	Rev_Size_yadA_Del	CGATTCCAAATATCTCTCGCAGG
	fliC deletion in I1/I2 and verification
	Fw_fliC	GGCATGATTATCCGTTTCTGC
	Fw_CmFRT_ForFliC	TTGGCGTTGCCGTCAGTCTCAGTTAATCAGGTTACAACGAATTTAAATGGCGCGCCTTAC
	Rev_CmFRT_ForFliC	CCCAATACGTAATCAACGACTTGCAATATAGGATAACGAATCGCCTACCTGTGACGGAAG
	Rev_fliC	TCCCAGCGATGAAATACTTGC

Table 4 .

 4 Genome assembly characteristics and metrics. Table showing data related to sequencing data, genome assembly quality and description.

	Case N°	Strain	Read number	Genome coverage average (SD)	Genome size (base pairs)	Contig number	GenBank Assembly accession	N50	L50	GC%	Complete BUSCOs	CDS	Potential plasmid gain	Potential plasmid loss	CDS gained	CDS lost	Mutation per base per year
	1	I1 R1	5,761,693 7,181,866	336 (170) 418 (200)	5,124,061 5,133,911	105 97	JARFPI000000000 JARFPF000000000	193,231 193,231	9 10	50.81 50.81	440/440 440/440	4,740 4,754	1	0	n.a 35	19	8.38x10 -6
	2	I2 R2	5,591,377 5,271,593	316 (144) 308 (132)	5,300,036 5,122,701	130 111	JARFPH000000000 JARFPE000000000	194,756 208,913	9 8	50.69 50.77	440/440 440/440	4,914 4,775	0	1	n.a 31	170	1.11x10 -6
	3	I3 R3	5,594,583 4,748,473	341 (138) 289 (118)	4,900,901 4,904,465	99 98	JARFPG000000000 JARFPD000000000	137,146 137,146	11 11	50.79 50.79	440/440 440/440	4,511 4,518	1	0	n.a 14	9	4.51x10 -6

Analysis of in-patient evolution of Escherichia coli reveals potential links to relapse of bone and joint infections

  Stanislas THIRIET-RUPERT 1 , Jérôme JOSSE 2 , David PEREZ-PASCUAL 1 , Jason TASSE 2 , Camille ANDRE 2 , Lélia ABAD 2 , David LEBEAUX 1 , Jean-Marc GHIGO 1 , Frédéric LAURENT 2 and Christophe BELOIN* 1

to 8.38x10 -6 ) as compared to what was reported in the literature (6.9x10 -7 for E. coli ED1a in

ACKNOWLEDGEMENTS.

We thank Dr. Emma Colucci-Guyon and Laurent Boucontet for critical reading of the manuscript. We are grateful to the Zebrafish Projects Hub of the Institut Pasteur and especially to Dr. Emma Colucci-Guyon and Laurent Boucontet who carried out the infection experiments in zebrafish larvae. We are grateful to Dr. Monica Rolando for kindly providing the antiflagellin antiserum. We are grateful to Emmanuel Chanard from Cerballiance for providing strains and microbiological information.

FUNDINGS

This work was supported by the French National Research Agency (ANR), project EvolTolAB (ANR-18-CE13-0010), by the French government's Investissement d'Avenir Program, Laboratoire d'Excellence "Integrative Biology of Emerging Infectious Diseases" (grant n°ANR-10-LABX-62-IBEID) and by the Fondation pour la Recherche Médicale (grant DEQ20180339185). S.T.-R was supported by the French National Research Agency (ANR), project EvolTolAB (ANR-18-CE13-0010).

DATA AVAILABILITY

All sequencing reads as well as genome assembly for each strain were deposited in NCBI under the BioProject accession number PRJNA932763. The code used in this study is available at https://github.com/Sthiriet-rupert/BJIs_adaptation.

COMPETING FINANCIAL INTERESTS

The authors declare no competing financial interests.

control, MG1655 strain (negative control) and CFT073 (positive control). Each image is representative of 5 biological replicates. Fifteen zebrafish were used for each experiment in C and D.

Supplementary Figures legends

Supplementary Figure 1. Optimization of competition protocols. (A) Assessment of the fitness cost associated to the expression of each fluorescence marker (GFP or mars) in each strain (n=3 biological replicates). There is a slight cost to produce mars marker in I3 and R3 strains as compared to the production of GFP. Therefore, both combinations of fluorescence marker were used in the competition assays for all strains. Statistical significance according to two-tailed Mann-Whitney test comparing initial and relapsing strains (*: p-value < 0.05, **: p-value < 0.01). (B) Comparison of FACS approach and Sanger based QSVanalyzer approach to assess the proprotion of each strain in a set of mixed samples with different proportions of I1 and R1 strains taged with GFP and mars, respectively (n=3 biological replicates). Supplemental Figure 2. Motility of each strain assessed on low agar plates at 30 and 37°C. A E. coli MG1655 strain deleted for fliE to fliR genes is used as a negative control and a E. coli MG1655 strain in which the motility is increased by the insertion of an IS1 in yeaJ/dgcJ encoding a diguanylate cyclase resulting in an increase of cdi-GMP, which in turn fuels motility was used as a positive control. Supplementary Figure 3. Influence of temperature on I1 and R1 strain motility. Motility of each strain at 30°C (A) or 37°C (B) when pre-cultured in liquid at 30°C or 37°C. (C) Motility of each strain deleted for fliC when precultured at 30°C or 37°C and incubated at 30°C for motility testing. (D) Production of flagellin (top panel) and RNA polymerase alpha subunit used as a loading control (bottom panel) in I1 and R1 strains assessed by Western blot analysis in different combinations of temperature in liquid cultures. Supplementary Figure 4. A. Auto-aggregation kinetic of all clinical strains. In each case, the infection strain is compared to the relapse strain as well as a strain of E. coli MG1655 in which flu is under the control of a constitutive promoter (PCL-flu) used as a positive control and a strain of E. coli MG1655 in which flu is deleted strain was used as a negative control and CFT073 strain as a positive control. Fifteen zebrafish were used for each experiment and no significant statistical different was identified between the infection strains and their relative relapse strains according to log-rank test for the survival and two-tailed Mann-Whitney test for the bacterial burden comparing initial and relapsing strains. Supplementary Figure 6. Growth capacities in different conditions. (A) Growth curve and related parameters inferred by the R package GrowthCurver (values for which the difference between both strains was significant (two-tailed Mann-Whitney test, p-value < 0.05) are highlighted in bold red) for I1 and R1 strains grown in M63B1 minimum medium supplemented with 0.4% pyruvate (n=5 biological replicates). (B-F) Growth curve and related parameters inferred by the R package GrowthCurver (values for which the difference between both strains was significant (two-tailed Mann-Whitney test, p-value < 0.05) are highlighted in bold red) for I2 and R2 strains grown in M63B1 minimum medium supplemented with 0.4% lactate (B), Maltose (C), Glucose (D) as well as in LB medium (E) and LB medium supplemented with 500 µM of 2,2'-Bipyridyl (F) (n=5 biological replicates).

Supplemental

Supplementary Results

Auto-aggregation and RDAR morphotypes on congo-red plates

Once settled and in order to colonize host environment, E. coli has a very large arsenal of adhesion factors. Some of them are known to promote biofilm formation via homotypic selfinteractions also called auto-aggregation [1,2], which could help facing stressful conditions such as the one encountered in vivo. The production of other factors such as curli and cellulose are also of critical importance [3,4] and induces a red, dry, and rough (rdar) morphotype when the strain is plated on LB agar supplemented with congo red. However, none of these phenotypes showed any difference when comparing the relapse strains to their related ancestors, all strains even having very weak auto-aggregation capacities (supplementary Figure S1.A.B).

Still, inter-couple differences were identified on congo red plates: (i) the first couple had a pink, dry and rough (pdar) morphotype characterizing the production of curli only; (ii) the second couple had a smooth and white morphotype (saw) showing the absence of matrix component production; and (iii) the third couple had pink and smooth morphotype (pas) showing the production of cellulose only.

Temperature specific motility

Because strain motility could have an impact on the virulence and biofilm forming capacity, we tested each strain on LB plates containing low amount of agar at 30°C and 37°C.

Motility is sometime more visible at lower temperature. The second and third couple strains had no motility (supplemental Figure S2). However, an interesting phenotype was observed for the first couple. Although both strains were motile, we observed some flares characteristic of the apparition of a hypermotile clone within the population and spreading faster than the others (supplementary Figure S2). Because this phenomenon was only observed at 37°C, we investigated the importance of temperature by testing all combinations between the temperature used to produce the liquid inoculum and the temperature the motility plates were incubated at.

We showed that regardless of the temperature used for the liquid culture, incubating the motility plate at 30°C always showed a homogenous motile phenotype, while plates inoculated with the same samples but incubated at 37°C always showed the apparition of flares (supplementary Figure S3.A.B.C). This suggests that the process is reversible and a westernblot analysis using anti-flagellin antibodies showed that both strains only produce flagellin at 30°C, even when pre-incubated at 37°C (supplementary Figure S3.D). Therefore, it seems that some mutants arise on plates incubated at 37°C in which the capacity to swim is restored.

Because the apparition of hyper-motile mutants suggests a high mutation rate, we used a Luria-Delbruck fluctuation test to evaluate it. Resistance to Rifampicin was used in these assays and strains I1 and R1 had a mutation rate of 1.58x10 -6 and 1.39x10 -6 mutations per generation, respectively, which was by far higher than MG1655 as estimated in the same assay (3.21x10 -8 mutations per generation), which was in accordance with what was reported in the literature for the same kind of experiment in Escherichia coli K12 MG1655 (see Figure 4 from [5]). Such a high mutation rate could be the result of a loss of function in genes responsible for DNA repair such as the mut genes. However, no such loss of function could be identified in the genome of the clinical strains.

Genome analysis

To confirm that in each case the relapse strain was the same as the initial infection one all of them were subjected to whole genome sequencing (WGS) and the genomic characteristics analyzed. Each case showed that both strains were indeed the same according to their phylogroup, serotype, ANI, MLST, FimH and FumC type (Table 1).

If the genome assemblies and completeness were good for all strain, the first notable difference identified between the infection strains and their related relapse strains was the number of gained and lost genes (supplementary Table S4). Some of these genes corresponded to a whole contig that matched known plasmids, therefore constituting potential gain or loss of plasmids (supplementary Data 2). R1 relapse strain gained a potential plasmid bearing Colicin E1 production genes, R2 relapse strain lost a potential plasmid containing hypothetical protein coding genes and R3 relapse strain gained a potential plasmid containing a Toxin/Antitoxin (T/A) system (HigA/HigB). However, PacBio sequencing technology would be required to rule out what proportion of the gene gains and losses was plasmid-mediated or the consequence of big chromosomic insertions or deletions.


Axel Osses
email: axosses@dim.uchile.cl

Nicolas Barnafi
email: nbarnafi@cmm.uchile.cl.

Nicol Ás

A Barnafi

FRANK-OSEEN TYPE SOLUTIONS WITH ORBITS AND VORTICES

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The nematic liquid crystal theory studies the structure of the principal directions of the molecules of a liquid crystal, which is an intermediate state between a solid and a liquid [START_REF] De Gennes | The physics of liquid crystals[END_REF]. Although this theory is typically used for studying electronic devices [START_REF] Bahadur | Liquid crystal-applications and uses[END_REF], other unconventional applications of ths theory have been established, such as proto-planet formation [START_REF] Rica | Formation de structures via une instabilité gravitationnelle[END_REF] and more recently the orientation of cardiac fibers, both experimentally [START_REF] Auriau | The nematic chiral liquid crystal structure of the cardiac myoarchitecture: disclinations and topological singularities[END_REF] and mathematically [START_REF] Barnafi | Modeling of cardiac fibers as oriented liquid crystals[END_REF]. Out work was motivated by the discovery of higher energy solutions of the Frank-Oseen equations on rings, computed in [START_REF] Barnafi | Modeling of cardiac fibers as oriented liquid crystals[END_REF]Sec. 2.2.3].

The simplest case of the Frank-Oseen equations of nematic liquid crystals [START_REF] Frank | On the theory of liquid crystals[END_REF] consists in finding a spatial distribution of unit vectors that are critical points of a certain global energy. The following is the simplest case of the Frank-Oseen theory under the one-constant hypothesis of lyquid crystal theory [START_REF] Majumdar | Landau-De Gennes theory of nematic liquid crystals: the Oseen-Frank limit and beyond[END_REF].

Definition 1. Given a regular domain Ω in R N , N = 2, 3, we say that d is a 2-nematic field if it is a critical point of the functional [START_REF] Armitage | Elliptic functions[END_REF] min

s.t.|d|=1 1 2 D |∇d| 2 - D f • d,
where f ∈ L 2 (Ω) is a given force and the solution is subject to the unit modulus constraint. That is, d is a unitary field solution of the following non linear system

∆d + 2λ 2 (d)d = -f, d • d = 1, in D, λ 2 (d) := 1 2 |∇d| 2 -f • d , (2)
where λ 2 the point-wise energy potential. If f = 0 we say that the 2-nematic field is homogeneous, non-homogeneous otherwise and conservative if f • d = 0 with total energy E 2 (Ω) = Ω λ 2 in the latter case. If the solution of (2) is in H 1 (Ω) we say this is a finite energy solution and an infinite energy solution otherwise.

The subscript 2 makes emphasis on the presence of the L 2 -norm for the gradient in (1), and will be instrumental for the definition of p-nematic liquid crystals in Section 5. In the classical theory [START_REF] Brezis | The interplay between analysis and topology in some nonlinear PDE problems[END_REF], we can consider f ∈ L 2 (Ω) N and seek for solution of (1) or [START_REF] Auriau | The nematic chiral liquid crystal structure of the cardiac myoarchitecture: disclinations and topological singularities[END_REF] with finite energy d ∈ H 1 (D) N , but we will build solutions that do not belongs to H 1 (D) N in a neighborhood of the origin.

There are other more complete theories for nematic crystals such as the Ericksen [START_REF] Ericksen | Liquid crystals with variable degree of orientation[END_REF], the Landau-De Gennes theories [START_REF] De Gennes | The physics of liquid crystals[END_REF][START_REF] Ball | Orientable and non-orientable line field models for uniaxial nematic liquid crystals[END_REF][START_REF] Ball | Nematic liquid crystals : from Maier-Saupe to a continuum theory[END_REF], and the Ginzburg-Landau theory of superconductors [START_REF] Ginzburg | On superconductivity and superfluidity (what i have and have not managed to do), as well as on the 'physical minimum' at the beginning of the 21 st century[END_REF][START_REF] Bethuel | Asymptotics for the minimization of a Ginzburg-Landau functional[END_REF].

Some of the solutions encountered in this paper were motivated from the geometry of the formation of planets and the distribution of their orbits, not only in our solar system but also in other planetary systems, geometries that exhibit orbits and vortices [START_REF] Bae | Structured distributions of gas and solids in protoplanetary disks[END_REF][START_REF] Dong | Multiple disk gaps and rings generated by a single super-Earth. ii. spacings, depths, and number of gaps, with application to real systems[END_REF]. This is part of the branch of Astronomy and Astrophysics that studies the dynamics of protoplanetary disks. We can think of a flat disk containing a fluid that rotates with an angular velocity and that under the action of gravity and magnetic fields is dragging matter until it accumulates in certain orbits. Even if the rigorous connection with minimization of energies of type [START_REF] Armitage | Elliptic functions[END_REF] or equations of type [START_REF] Auriau | The nematic chiral liquid crystal structure of the cardiac myoarchitecture: disclinations and topological singularities[END_REF] is not yet completely clear, there are some similarities with solutions obtained by perturbation analysis of a non-linear fluid accretion model in [START_REF] Rica | Formation de structures via une instabilité gravitationnelle[END_REF].

The shape and structure of these disks where planets are forming can be observed and studied with radio telescopes. Numerical simulations can also be performed from the equations of magneto-hydrodynamics [START_REF] Bae | Structured distributions of gas and solids in protoplanetary disks[END_REF] where after multiple spins we arrive to structures similar to those observed in radio-astronomy. The problem of fit of planet orbits, autosimilar and quantified solutions appears also in the context of the scale relativity theory [START_REF] Nottale | Scale relativity and fractal space-time. A new approach to unifying relativity and quantum mechanics[END_REF].

There is also a recent connection between nematic theory and the distribution of the orientation of muscular cardiac fibers, both from the theoretical, numerical [START_REF] Barnafi | Modeling of cardiac fibers as oriented liquid crystals[END_REF] and the experimental point of view [START_REF] Auriau | The nematic chiral liquid crystal structure of the cardiac myoarchitecture: disclinations and topological singularities[END_REF], where Frank-Oseen models naturally appears in the presence of external forces in the case of a single ventricle [START_REF] Barnafi | Modeling of cardiac fibers as oriented liquid crystals[END_REF]. Again, the rigorous link with minimization of energies of type [START_REF] Armitage | Elliptic functions[END_REF] is not yet completely clear.

Motivated from geometries from these applied fields, we build novel solutions of (2) that have infinite energy [START_REF] Brezis | Harmonic maps with defects[END_REF]. The solutions are exact in two dimensions (i.e. with f = 0) and approximate in the three dimensional case (i.e. with f = 0 but with a certain radial power decay from the origin). We focus only on the cases where the domain D is either a disk, a cylinder, and a sphere.

The rest of the paper is organized as follows. In Section 2 we compute analytic solutions on a disk. In Section 3 we extend the disk solutions to a cylinder. This formulation gives rise to a non-homogeneous problem, so we show three types of solutions where we try to maximize the power decay of the external force. In Section 4 we perform a similar study but for the case of a spherical domain. In Section 5 we define p-nematic liquid crystals for p ≤ 2 and establish connections with the standard p = 2 model. Finally, we discuss our results in Section 6.

Exact solutions in a two-dimensional disk

Let us consider a disk of radius R centered at the origin of coordinates which we will represent in polar coordinates (r, θ) given by D = {|r| ≤ R}. Problem (2) can be rewritten as the solution of the following nonlinear Frank-Oseen equation: find a unit vector field d such that:

(3) ∆d + λd = d rr + 1 r d r + 1 r 2 d θθ + 2λ 2 d = 0 in D, where (4) λ 2 = 1 2 |∇d| 2 = 1 2 |d r | 2 + 1 r 2 |d θ | 2 .
We will look for solutions of this equation with boundary conditions on the boundary ∂D of D, yet to be defined. As we show later, we will not worry about possible singularities at r = 0-that is, we can obtain infinite energy solutions-because we can always remove a small neighborhood of the origin. Henceforth it will be useful for us to use the rotation matrix at an angle α and two related matrices:

(5)

Q α = cos α sin α -sin α cos α , I = Q 0 = 1 0 0 1 , I ⊥ = Q π/2 = 0 1 -1 0 ,
which satisfy the following properties: dQ dα = I ⊥ Q, and I ⊥ I ⊥ = -I. These properties will be used extensively throughout this work.

In two dimensions, we will look for a solution with a rotating boundary condition

(6) e(nθ) = -sin nθ cos nθ if r = R,
so we propose the following solution ansatz :

(7) d(r, θ) = Q α(r) e(nθ),
where α(r) is a function to be determined. Leveraging the properties of the rotation matrix, we readily see that

d r = α I ⊥ d, d rr = (-|α | 2 I + α I ⊥ )d, d θ = -nI ⊥ d, d θθ = -n 2 d, (8)
which in (3) yields ( 9)

-|α | 2 I + α I ⊥ + α r I ⊥ - n 2 r 2 I + 2λI = 0, λ 2 = 1 2 |α | 2 + n 2 r 2 .
We highlight that the components of I cancel out. By canceling also the components of I ⊥ and using the boundary condition we obtain

   α + α r = 0 r ∈ (0, R), α(R) = 0,
from where we obtain the general solution for α:

(10) α = c 1 log(r/R), λ 2 = c 2 1 + n 2 2r 2 ,
where c 1 is an arbitrary constant. Finally, we established the following result.

Theorem 2. For each n ≥ 1, there exist a family of solutions of the Frank-Oseen equation:

(11)      ∆d + |∇d| 2 d = 0 in D, d • d = 1 in D, d = e(nθ) on ∂D.
in the disk D = {(r, θ), |r| ≤ R} where e was defined in [START_REF] Ball | Orientable and non-orientable line field models for uniaxial nematic liquid crystals[END_REF], given by the unitary field:

(12) d(r, θ) = Q c 1 log(r/R) e(nθ)
and where c 1 ∈ R is an arbitrary constant and Q α is a rotation matrix defined in [START_REF] Ball | Nematic liquid crystals : from Maier-Saupe to a continuum theory[END_REF].

See Figures 1 and2 for a graphical representation of the family of solutions (12) of Theorem 2. We note that the solutions exhibits for n = 1 alternate attractive and repulsive orbits in rings with exponentially growing separation, for n ≥ 2 the attractors are spirals. Remark 3. It is easy to verify that the solutions (12) are self-similar with respect to the origin. Indeed, if d(r, θ) is the solution given by (12) then d k (r , θ ) = d(r, θ) for r = r/a k , θ = θ and a k = R exp 2kπ c 1 are also solutions of the same equations (2) with the same boundary conditions for r = 1. In this ring, it is not difficult to check that the energy of the solution using (10) is given by

E 2 (D \ D ε ) = 2π 0 R ε λ 2 (r, θ)r drdθ = π(c 2 1 + n 2 ) log(R/ε),
so we see that the H 1 -norm of the solution tends to infinity as ε → 0. We note that if we fix n = 1 and c 1 = 0, we recover a minimal energy solution. For other values of c 1 and n, the solutions are critical points of the energy functional (1) (with f = 0) with higher energy levels. See Section 5 for a more detailed discussion about the energy of the solutions.

Approximate solutions in a cylinder

In this section, we extend the solutions obtained in a disk to the case of a cylinder. For this we will consider the cylinder in cylindrical coordinates:

C = {(r, θ, z) | 0 ≤ r ≤ R, 0 ≤ θ < 2π, 0 ≤ z ≤ H}.
The Frank-Oseen equation in polar coordinates is written as searching for a unitary field d such that:

(13) d rr + 1 r d r + 1 r 2 d θθ + d zz + 2λ 2 d = -f in C, with some right hand side f such that f • d = 0, where (14)
λ 2 = 1 2 |d r | 2 + 1 r 2 |d θ | 2 + |d z | 2 .
We define Notice that e • e ⊥ = 0 and | e| = 1, | e ⊥ | = 1. We also introduce the rotation matrices in the xy plane of the form ( 16)

Q α =   cos α sin α 0 -sin α cos α 0 0 0 1   , I 12 =   1 0 0 0 1 0 0 0 0   , I ⊥ 12 =   0 1 0 -1 0 0 0 0 0   , I 3 =   0 0 0 0 0 0 0 0 1   ,
and I the usual identity matrix in R 3 . As before, we will make use of the following properties:

Q α Q t α = I, I 12 I ⊥ 12 = 0, I ⊥ 12 I ⊥ 12 = -I 12 .
Additionally, Q α and I ⊥ 12 conmute. For the cylinder case, we consider a solution ansatz of the form

d(r, θ, z) = Q α(r) ê(nθ, β(ξ)),
where we will further consider three types of solution. One with a height rotation ξ = z and two with an elevation rotation ξ = z/r. The difference between the last two will reside in the treatment of the forcing term.

3.1. Cylindrical case 1. In this subsection, we look for solutions of the form

d(r, θ, z) = Q α(r) e (nθ, β(z)) ,
where α(r) and β(z) are functions to determine, with α(R) = 0 and β(0) = 0 as suggested by the boundary conditions. The partial derivatives are given as follows:

d r = α I ⊥ 12 d, d rr = (-|α | 2 I 12 + α I ⊥ 12 )d, d θ = -nI ⊥ 12 d, d θθ = -n 2 I 12 d, d z = β Q α e ⊥ , d zz = -|β | 2 d + β Q α e ⊥ .
By replacing them in [START_REF] De Gennes | The physics of liquid crystals[END_REF] we obtain ( 17)

-|α | 2 + n 2 r 2 I 12 d -|β | 2 d + α + α r I ⊥ 12 d + β Q α e ⊥ + 2λ 2 d = -f 1 , λ 2 = 1 2 cos 2 β |α | 2 + n 2 r 2 + |β | 2 , f 1 = g 1 Q α e ⊥ , g 1 (r, z) = |α | 2 + n 2 r 2 sin β cos β,
where f 1 comes from the terms that do not vanish on the left hand side and can be rewritten as:

f 1 = |α | 2 + n 2 r 2 (sin 2 βI 12 -cos 2 βI 3 )d = sin β cos β |α | 2 + n 2 r 2 Q α e ⊥ .
Notice that d, I ⊥ 12 d and Q α e ⊥ are mutually orthogonal since

I ⊥ 12 d • d = 0, Q α e ⊥ • d = Q α e ⊥ • Q α e = e t Q t α Q α e ⊥ = e • e ⊥ = 0,
and

I ⊥ 12 d • Q α e ⊥ = I ⊥ 12 Q α e • Q α e ⊥ = Q t α I ⊥ 12 Q α e • e ⊥ = Q t α Q α I ⊥ 12 e • e ⊥ = I ⊥ 12 e • e ⊥ = 0.
This implies in particular that f 1 • d = 0. Equating the components of I ⊥ 12 and Q α e ⊥ on the left hand side to zero, we obtain (18)

           rα + α = 0 r ∈ (0, R), β = 0 z ∈ (0, Z), β(0) = 0, α(R) = 0, from which it follows that (19) α = c 1 log(r/R), β(z) = c 2 z,
for arbitrary real constants c 1 and c 2 . We have thus obtained the following result.

Theorem 5. For each n ≥ 1, family of solutions of the non-homogeneous Frank-Oseen equations

(20) ∆d + |∇d| 2 d = -f 1 in C, d • d = 1 in C, in the cylinder C = {(r, θ, z) | 0 ≤ r ≤ R, 0 ≤ θ < 2π, 0 ≤ z ≤ H} with right hand side f 1 such that f 1 • d = 0 and f 1 of order O(1/r 2 )
, is given by

(21) d(r, θ, z) = Q c 1 log(r/R) e (nθ, c 2 z)
where c 1 , c 2 ∈ R are arbitrary constants, Q α is the rotation matrix defined in ( 16) and e was defined in (15).

Cylindrical case 2.

In the previous calculations, it would be interesting to absorb f 1 by solving β + g 1 = 0, but the solution of this equation is impossible since β only depends on z and g 1 depends on (r, z). This suggest to take β = β(z, r). We will try the case β = β(z/r), i.e. function of the elevation angle-but other choices could be possible-, that is, search now for a solution of type

(22) d(r, θ, z) = Q α(r) e (nθ, β(z/r)) ,
where α(r) and β(ξ), ξ = z/r such that α(R) = 0 and β(0) = 0 are functions to determine. The partial derivatives are now given by

d r = α I ⊥ 12 d - z r 2 β Q α e ⊥ , d rr = (-|α | 2 I 12 + α I ⊥ 12 )d + z 2 r 4 β + 2 z r 3 β Q α e ⊥ - z 2 r 4 |β | 2 d - 2z r 2 α β I ⊥ 12 Q α e ⊥ , d θ = -nI ⊥ 12 d, d θθ = -n 2 I 12 d, d z = 1 r β Q α e ⊥ , d zz = - 1 r 2 |β | 2 d + 1 r 2 β Q α e ⊥ .
By replacing them in [START_REF] De Gennes | The physics of liquid crystals[END_REF] we obtain (23)

-|α | 2 + n 2 r 2 I 12 d - 1 r 2 1 + z 2 r 2 |β | 2 d + α + α r I ⊥ 12 d + 1 r 2 1 + z 2 r 2 β + z r β Q α e ⊥ + 2λ 2 d = -f 1 -f 2 λ 2 = 1 2 cos 2 β |α | 2 + n 2 r 2 + 1 r 2 1 + z 2 r 2 |β | 2 f 1 = g 1 Q α e ⊥ , g 1 (r, z) = |α | 2 + n 2 r 2 sin β cos β f 2 = - 2z r 2 α β I ⊥ 12 Q α e ⊥
and the two terms f 1 and f 2 on the right hand satisfy

f 1 • d = f 2 • d = 0.
The last equality comes from the fact that

I ⊥ 12 Q α e ⊥ • d = Q α I ⊥ 12 e ⊥ • Q α e = Q t α Q α I ⊥ 12 e ⊥ • Q α e = I ⊥ 12 e ⊥ • e = 0.
Equating the components of I ⊥ 12 and of Q α e ⊥ on the left hand side to zero, we obtain (24)

           rα + α = 0 r ∈ (0, R), (1 + ξ 2 )β + ξβ = 0 ξ ∈ (0, ∞), α(R) = 0 , β(0) = 0 .
Note that the domain of ξ was computed from ξ = z/r as z → 0 and as r → 0. The solution of ( 24) is given by

(25) α = c 1 log(r/R), β = c 2 sinh -1 (z/r),
where sinh -1 (ξ) = log ξ + 1 + ξ 2 . Thus, we have obtained the following result.

Theorem 6. For each n ≥ 1, there is a family of solutions of the non-homogeneus Frank-Oseen equation

(26) ∆d + |∇d| 2 d = -f 1 -f 2 in C, d • d = 1 in C, in the cylinder C = {(r, θ, z) | 0 ≤ r ≤ R, 0 ≤ θ < 2π, 0 ≤ z ≤ H}, with right hand side such that f 1 • d = 0, f 2 • d = 0 and f 1 = O(1/r 2 ) and f 2 = O(z/r 3 ), is given by (27) d(r, θ, z) = Q c 1 log(r/R) e nθ, c 2 sinh -1 (z/r) ,
where c 1 , c 2 ∈ R are arbitrary constants, Q α is the rotation matrix defined in ( 16) and e was defined in (15).

3.3. Cylindrical case 3. We can go further and search for a solution that absorbs f 1 (but not f 2 ). Equating the components of I ⊥ 12 on the left hand side to zero and the components of Q α e ⊥ on the left hand side to g 1 , we obtain

           rα + α = 0 r ∈ (0, R), (1 + ξ 2 )β + ξβ + r 2 |α | 2 + n 2 sin β cos β = 0 ξ ∈ (0, ∞) α(R) = 0 , β(0) = 0 .
The solution of the first equation is again given by α(r) = c 1 log(r/R). Therefore,

(1 + ξ 2 )β + ξβ + A 2 sin β cos β = 0, β(0) = 0, where A 2 = c 2 1 + n 2 . Now making the change of variables γ(v) = β(sinh v), we see that γ solves γ + A 2 sin γ cos γ = 0, γ(0) = 0, whose solution is γ(v) = am(c 2 v; k 2 )
, where am is the Jacobi elliptic amplitude function with parameter (modulus) 0 < k 2 < 1 and c 2 = k/A [START_REF] Armitage | Elliptic functions[END_REF]. This is because the am function satisfies the second order equation am (ϕ) + k 2 sn(ϕ) cn(ϕ) = 0, sin(am(ϕ)) = sn(ϕ) and cos(am(ϕ)) = cn(ϕ). From this it follows that α = c 1 log(r/R), (28)

β = am(ϕ), ϕ = c 2 sinh -1 (z/r), c 2 = k c 2 1 + n 2 . ( 29
)
If we use that am (ϕ) = dn(ϕ) then

f 2 = - 2c 1 z r 3 dn(ϕ)Q α I ⊥ 12 e ⊥ (nθ, am(ϕ))
that satisfies f 2 • d = 0 and because cn, sn and dn are bounded, we have that f 2 is of order O(z/r 3 ). Thus we have obtained the following result.

Theorem 7. For each n ≥ 1, there is a family of solutions of the non-homogeneus Frank-Oseen equation

(30) ∆d + |∇d| 2 d = -f 2 in C, d • d = 1 in C, in the cylinder C = {(r, θ, z) | 0 ≤ r ≤ R, 0 ≤ θ < 2π, 0 ≤ z ≤ H}, with right hand side f 2 such that f 2 • d = 0 and f 2 = O(z/r 3 ) given by (31) d(r, θ, z) = Q c 1 log(r/R) e nθ, am(c 2 sinh -1 (z/r); k 2 ) , c 2 = k c 2 1 + n 2
, where c 1 , c 2 ∈ R are arbitrary constants, Q α is the rotation matrix defined in (16), am(•; k 2 ) is the Jacobi amplitude of modulus k 2 ∈ (0, 1) and e was defined in [START_REF] Ericksen | Liquid crystals with variable degree of orientation[END_REF].

Remark 8. The solutions in (31) can be easily computed using the Jacobi elliptic functions sn and cn because e defined in (15) satisfies Remark 9. It is natural to ask why we have not provided the computations for the fully harmonic case in which both load vectors f 1 and f 2 are absorbed through the scalar functions as in the previous cases. In fact, this is not possible, because the differential equations obtained in such case are given by

r 2 α + rα + 2ξrα β tan β = 0, α(R) = 0 (1 + ξ 2 )β + ξβ + r 2 |α | 2 + n 2 sin β cos β = 0, β(0) = 0, with ξ = z/r. ( 33
)
Performing the changes of variables

δ(u) = α(e u ), u = log(r), γ(v) = β(sinh v), v = sinh -1 (ξ)
we obtain the following system of second order non linear ordinary differential equations

δ + 2δ γ tan γ = 0, δ(log R) = 0, γ + |δ | 2 + n 2 sin γ cos γ = 0, γ(0) = 0. ( 34
)
We note that this implies that α the variables depend on the argument of the other variable, which contradicts the form of the ansatz considered. Indeed, previous computations were successful because they allowed for a decoupling of the scalar functions α and β.

Remark 10. Notice that the singularity at the origin in all the cases of this section is of type log(r) as for the two dimensional case, but around all the vertical axis z = 0, so it is a vortex [START_REF] Sandier | Vortices in the magnetic Ginzburg-Landau Mmodel[END_REF]. We could also remove this singularity by considering the equations in a perforated cylinder

C \ C ε with C = {(r, θ, z) | |r| < ε, 0 ≤ θ < 2π, |z| ≤ H}.
Otherwise, we will see in Section 5 that these solutions are not in H 1 (C) but in W 1,p (C) for p < 2.

We have numerically computed the solutions of equations ( 20), (26), and (30), considering a symmetric cylinder with z ∈ (-H, H). The 3D plot of the streamlines obtained from the solution is shown in Figure 3, where we highlight the presence of the orbits (n = 1) or spirals (n ≥ 2) in the xy-plane and a vortex in the z axis for all the solutions. Instead, in Figure 5 we have plotted the first component of the solution in three different settings: the mid-xy plane, its sum on the same plane and then the solution on the mid-xz plane. Interestingly, self-similar structures are observable aloso along the xz plane (see the third column of Figure 5). In Figure 4 we provide the observation of the structure of a proto-planetary disk to see the similarities that motivated us to build this type of solutions.

Approximate solutions in a ball

Let us now consider the ball in spherical coordinates

B = (r, θ, φ) | |r| ≤ R, 0 ≤ θ < 2π, - π 2 ≤ φ < π 2
where θ is the azimuth angle and φ is the elevation angle, i.e., φ = 0 corresponds to the mid horizontal plane. The non-homogeneous Frank-Oseen equation in spherical coordinates given by searching for a unitary field d such that:

(35) d rr + 2 r d r + 1 r 2 cos 2 φ d θθ + 1 r 2 cos φ (cos φ d φ ) φ + 2λ 2 d = -f in B
where f • d = 0 and (36)

λ 2 = 1 2 |d r | 2 + 1 r 2 cos 2 φ |d θ | 2 + 1 r 2 |d φ | 2 .
As in the cylinder case, we separate our cases according to the treatment of the right hand side. where α(r) and β(φ) are functions to be determined such that α(R) = 0, β(0) = 0, that is, the vector field rotates horizontally in the central plane xy. The partial derivatives are

d φ = β Q α e ⊥ , d φφ = -|β | 2 d + β Q α e ⊥ .
Replacing in (35) yields (38)

-|α | 2 I 12 d + α + 2α r I ⊥ 12 d - n 2 r 2 cos 2 φ I 12 d - |β | 2 r 2 d + 1 r 2 (β -tan(φ)β )Q α e ⊥ + 2λ 2 I = -f 3 -f 4 , λ 2 = 1 2 cos 2 β |α | 2 + n 2 r 2 cos 2 φ + 1 r 2 |β | 2 , f 3 = g 3 Q α e ⊥ , g 3 = sin β cos β|α | 2 , f 4 = g 4 Q α e ⊥ , g 4 = sin β cos β n 2 r 2 cos 2 φ .
Equating the components of I ⊥ 12 and of Q α e ⊥ on the left hand side, we obtain

             rα + 2α = 0 r ∈ (0, R), β -tan(φ)β = 0 φ ∈ - π 2 , π 2 , α(R) = 0, β(0) = 0,
from which we have that

α = c 1 1 R - 1 r , (39)
β = c 2 log 1 + sin φ cos φ , (40)
We have thus obtained the following result.

Theorem 11. For each n ≥ 1 there is a family of solutions of the non-homogeneous Frank-Oseen equations

(41) ∆d + |∇d| 2 d = -f 3 -f 4 in B, d • d = 1 in B, in the ball B = (r, θ, φ) | |r| ≤ R, 0 ≤ θ < 2π, -π 2 ≤ φ < π 2 with f 3 • d = 0, f 4 • d = 0, f 3 = O(1/r 2 ) and f 4 = O(1/r 4 ) of the form (42) d(r, θ, φ) = Q c 1( 1 R -1 r ) e nθ, c 2 log 1 + sin φ cos φ ,
where c 1 , c 2 ∈ R are arbitrary constants, Q α is the rotation matrix defined in ( 16) and e was defined in [START_REF] Ericksen | Liquid crystals with variable degree of orientation[END_REF].

Remark 12. For the solutions (42) of the spherical case in the previous theorem, they do not present orbits in the mid xy-plane, but only spirals.

Spherical case 2.

If we want to absorb f 3 from (35) we would need to solve (43)

             rα + 2α = 0 r ∈ (0, R), cos 2 φβ -sin φ cos φβ + n 2 sin β cos β = 0 φ ∈ - π 2 , π 2 , α(R) = 0 , β(0) = 0 .
The solution for α is the same as for the previous case. For β, if we define

(44) v = h(φ) = log 1 + sin φ cos φ with h(0) = 0, now if we make the change of variables γ(v) = β(h -1 (v)) then γ = cos 2 φβ -sin φ cos φβ and γ is a solution of (45) γ + n 2 sin γ cos γ = 0, γ(0) = 0.
Again, the solution is given by γ(v) = am(c 2 v; k 2 ) of modulus k, where c 2 = n/k. From the previous analysis we have that

α = c 1 1 R - 1 r , ( 46
)
β = am c 2 log 1 + sin φ cos φ ; k 2 , c 2 = k n , (47)
and we obtain the following result.

Theorem 13. For each n ≥ 1 there is a family of solutions of the non-homogeneous Frank-Oseen equations

(48) ∆d + |∇d| 2 d = -f 4 in B, d • d = 1 in B, in the ball B = (r, θ, φ) | |r| ≤ R, 0 ≤ θ < 2π, -π 2 ≤ φ < π 2 with f 4 • d = 0, f 4 = O(1/r 4 ) of the form (49) d(r, θ, φ) = Q c 1( 1 R -1 r ) e nθ, am c 2 log 1 + sin φ cos φ ; k 2 , c 2 = k n ,
where c 1 , c 2 ∈ R are arbitrary constants, Q α is the rotation matrix defined in (16), am(•; k 2 ) is the Jacobi amplitude of modulus k 2 ∈ (0, 1) and e was defined in [START_REF] Ericksen | Liquid crystals with variable degree of orientation[END_REF].

Remark 14. Again, all the solutions (49) of the spherical case in the previous theorem do not exhibit orbits in the mid xy-plane, but only spirals. See Figure 8 and compare with Figure 2.

Remark 15. Note that the singularity of the solution at the origin in the spherical cases are of type 1/r. There is also a singularity on the vertical axis (elevation φ = ±π/2) of type log π 2 -|φ| which is a vortex [START_REF] Sandier | Vortices in the magnetic Ginzburg-Landau Mmodel[END_REF]. This means the solutions are not in H 1 (B). We could remove these singularities on the origin and the vertical axis by considering the ball

B minus B ε = {(r, θ, φ) | |r| < ε, 0 ≤ θ < 2π, |φ| < arctan(R/ε)}.
Remark 16. All of the missing combinations of terms that can be absorbed into the left hand side yield the same contradiction detailed in Remark 9. We omit further details.

We have numerically computed the solutions of ( 41) and ( 48), with its streamlines shown in Figure 6, where the presence of a radial spiral instead of the logarithmic orbits can be appreciated. In Figure 7 we plot the first component of both solutions in the xy and xz planes. Interestingly, the solution quickly fades to the origin, which prevents the formation of orbits. It is also possible to see the discontinuity along the ±π/2 elevations. We conclude this by showing how the solution changes in the partially non-homogeneous case as the angle parameter n changes in Figure 8. We stretch the analogy of liquid crystals and protoplanets formation to show that, only anecdotically, the solutions resemble protoplanetary jet observations, as the one shown in Figure 4.

p-nematic fields

In the previously constructed solutions, we have obtained that they satisfy the Frank-Oseen equations ( 2), but they have infinite energy in H 1 (Ω), where Ω is either a disk in two dimension or a cylinder or a sphere in three dimensions. It is natural to ask if there is a modified potential which is (i) finite and (ii) such that the solutions are critical points (or minima). Given that the functions have been observed to belong to W 1,p (Ω), we define a vector field to be p-nematic by the following definition. Definition 17. Given a regular domain Ω ∈ R N , N = 2, 3, we say that d is a p-nematic liquid crystal in W 1,p (Ω), if given f in the dual space of W 1,p (Ω), it is a critical point of:

(50) min

s.t.|d|=1 1 p Ω |∇d| p - Ω f • d for p ≥ 1 in a certain domain Ω.
That is, a unit vector field d that satisfies the generalized Frank-Oseen equations:

∆ p d + pλ p (d)d = -f, d • d = 1 in Ω, λ p (d) := 1 p |∇d| p -|∇d| p-2 f • d , (51)
where

(52) ∆ p d = ∇ • (|∇d| p-2 ∇d)
is the p-Laplacian [START_REF] Evans | Partial Differential Equations[END_REF]. We say that the p-nematic field is homogeneous if f = 0, nonhomogeneous otherwise and conservative if f • d = 0 with total energy E p (Ω) = Ω λ p in the latter case.

Notice that the case p = 2 corresponds to the 2-nematic solutions encountered in the previous sections and introduced in Definition 1, but we recall that they had finite energy only away from the the origin or the vertical axis. We have the following purely algebraic relationship between the p = 2 case and the general case that will be useful in the following.

Lemma 18. Let d be a vector field in R N , then

(53) λ p = 2 p |∇d| p-2 λ 2 , which in the conservative case f • d = 0 reduces to (54) λ p = 2 p/2 p λ p/2 2 ,
and the corresponding energy in the conservative case is given by

(55) E p (Ω) = Ω λ p = 2 p/2 p Ω λ p/2 2 .
In the following, we provide estimates of the p-energy for most of the solutions obtained in the previous sections.

Theorem [START_REF] Majumdar | Landau-De Gennes theory of nematic liquid crystals: the Oseen-Frank limit and beyond[END_REF]. The energy of the solutions of the Frank-Oseen solutions of the previous sections are estimated by:

• Disk case: The solutions of Theorem 2 in the disk D of radius R > 0 do not belong to H 1 (D) and they belong to W 1,p (D) for all 1 ≤ p < 2. More precisely, if

D ε = {(r, θ), |r| < ε, 0 ≤ θ < 2π}, then E 2 (D) = +∞ E 2 (D \ D ε ) = π(c 2 1 + n 2 ) log(R/ε), E p (D) = 2π p(2 -p) (c 2 1 + n 2 ) p/2 R 2-p .
• Cylindrical cases: The solutions of Theorems 5, 6 and 7 in the cylinder C of radius R > 0 and height 2H > 0 do not belong to H 1 (C) and they belong to W 1,p (C) for all 1 ≤ p < 2. More precisely, if

C ε = {(r, θ, z), |r| < ε, 0 ≤ θ < 2π, |z| ≤ H}, then E 2 (C) = +∞ E 2 (C \ C ε ) = O (log(R/ε)) as ε → 0, E p (C) = O R 2-p 2 -p as p → 2.
• Spherical cases: The solutions of Theorems 11 and 13 in the ball B of radius R > 0 do not belong to H 1 (B) and they belong to W 1,p (B) for all 1 ≤ p < 3/2. More precisely, if

B ε,δ = {(r, θ, φ), |r| < ε, 0 ≤ θ < 2π, |φ| < π/2 -δ}, then E 2 (B) = +∞ E 2 (B \ B ε,δ ) = O 1 ε - 1 R log(1/δ) as ε → 0 and δ → 0, E p (B) = O R 3-2p 3 -2p as p → 3/2.
Proof. In the two dimensional case of a disk, it is easy to verify from ( 9), [START_REF] Brezis | Harmonic maps with defects[END_REF] and that

E 2 (D \ D ε ) = D\Dε λ 2 = π R ε c 2 1 + n 2 r 2 p/2 r dr = π(c 2 1 + n 2 ) log(R/ε)
where we see that as ε → 0 the integral diverges to +∞ and using the previous proposition for 1 ≤ p < 2

E p (D) = 2 p/2 p D λ p/2 2 = 2π p R 0 c 2 1 + n 2 r 2 p/2 r dr = 2π p(2 -p) (c 2 1 + n 2 ) p/2 R 2-p
from which we obtain the statement of the theorem for this case.

In the cylindrical case 1, using ( 17) and ( 19), we see that the energy in the cylinder C is given by

E p (C) = 2 p/2 p 2π 0 H -H R 0 λ p/2 2 r drdzdθ = 2π p H -H R 0 | cos(c 2 z)| p c 2 1 + n 2 r 2 + c 2 2 p/2 r drdz ≤ 4πH p (c 2 1 + n 2 + c 2 2 R 2 ) p/2 2 -p R p-2 .
For 1 ≤ p < 2 we have x p ≥ x 2 for x ∈ [0, 1] and H -H cos 2 (c 2 z)dz = H (1 + sinc(2c 2 H)), using also that |x| p/2 is a non decreasing function we have that

E p (C) ≥ 2πH(1 + sinc(2c 2 H)) p (c 2 1 + n 2 ) p/2 2 -p R p-2 .
In the same way, it is easy to check that

E 2 (C \ C ε ) ≤ 2πH(c 2 1 + n 2 + c 2 R 2 ) log(R/ε) and E 2 (C \ C ε ) ≥ πH(1 + sinc(2c 2 H))(c 2 1 + n 2 ) log(R/ε).
In the cylindical case 2, from ( 23) and ( 25), the energy in the cylinder C for 1 ≤ p < 2 is given by

E p (C) = 2π p H -H R 0 cos 2 (β(z/r)) c 2 1 + n 2 r 2 + |β (z/r)| 2 1 r 2 + z 2 r 4 p/2 rdrdz ≤ 4πH p R 0 c 2 1 + n 2 + 2c 2 2 r 2 p/2 rdr ≤ 4πH p(2 -p) (c 2 1 + n 2 + 2c 2 2 ) p/2 R 2-p ,
where β(ξ) = c 2 sinh -1 (ξ), and we used that | cos(β)| ≤ 1 and ξ(sinh -1 ) (ξ) = ξ(1 + ξ 2 ) -1/2 ≤ 1 for ξ = z/r. It is also easy to see that since

(1 + ξ 2 ) -1/2 ≥ (1 + (H/R) 2 ) -1/2 then E p (C) ≥ 4πH p(2 -p) (1 + (R/H) 2 ) -p/2 R 2-p .
In the same way, it is easy to check in this case that

E 2 (C \ C ε ) ≤ 2πH(c 2 1 + n 2 + 2c 2 ) log(R/ε) and E 2 (C \ C ε ) ≥ 2πH(1 + (R/H) 2
) -1 log(R/ε). For the cylindrical case 3, in which β is given this time by (28)), the analysis is analogous to the previous one by replacing β(ξ) = am(c 2 sinh -1 (ξ)), because in this case β (ξ) = c 2 dn(c 2 sinh -1 (ξ))(1 + ξ 2 ) -1/2 and the Jacobi function dn is strictly positive and more precisely has range in [ √ 1 -k 2 , 1]. We obtain thus the same upper bounds of the previous case and the same lower bounds up to the factor √ 1 -k 2 . For the spherical case 1, given by ( 38) and (39), the energy in the ball B (since in this case β (φ) = c 2 / cos φ) is given by

E p (B) = 2π p π/2 -π/2 R 0 c 2 1 cos 2 β(φ) r 4 + n 2 cos 2 β(φ) + c 2 2 r 2 cos 2 φ p/2
r 2 cos φ drdφ.

We can bound the energy from above by

E p (B) ≤ 2π p max{c p 1 , (n 2 + c p 2 ) p } π/2 -π/2 R 0 1 + R 2 r 4 cos 2 φ p/2 r 2 cos φ drdφ, ≤ 2π p max{c p 1 , (n 2 + c p 2 ) p }(1 + R 2 ) p/2 R 0 r 2-2p dr π/2 -π/2
(cos φ) 1-p dφ, which gives that the radial integral is finite for p < 3/2, and for the angular integral for the elevation angle φ we have

π/2 -π/2 (cos φ) 1-p dφ = 2 π/2 0 (sin ψ) 1-p dψ
and sin ψ ≈ ψ for ψ small so the integral is finite (say a constant C 3 > 0) for p < 2. Finally

E p (B) ≤ C 3 2π p max{c p 1 , (n 2 + c p 2 ) p }(1 + R 2 ) p/2 R 3-2p 3 -2p .
For a lower bound, we have that

E p (B) ≥ 2πc p 1 p R 0 r 2-2p dr π/2 -π/2 cos 2 (β(φ)) cos φ dφ,
and the integral in φ is strictly positive, since its argument is non negative, continuous and thus positive at some point φ = φ 0 with -π/2 < φ 0 < π/2, so we conclude by continuity that there exists a constant C 4 > 0 such that

E p (B) ≥ C 4 2πc p 1 p R 3-2p 3 -2p .
In the same way, using the same arguments as before, it is easy to check that there exist positive constants C 5 and C 6 such that

E 2 (B \ B ε,δ ) ≤ C 5 1 ε - 1 R log(1/δ) and E 2 (B \ B ε,δ ) ≥ C 6 1 ε - 1 R log(1/δ).
The factor log(1/δ) comes from the fact that

π/2-δ -π/2+δ (cos φ) -1 dφ = 2 π/2-δ 0 (sin ψ) -1 dψ = O (log(1/δ)) .
The proof for spherical case 2 given as in the precedent case by (38) but only changing the function β by (46), is analogous to the previous one.

One may naturally wonder if a 2-nematic solution behaves also as a p-nematic solution, under certain circumstances. We have the following.

Theorem 20. If d is a non-homogeneous 2-nematic solution associated to f with f • d = 0, then d is also a p-nematic non homogeneous solution associated to f defined by

(56) f = |∇d| p-2 f -(p -2)|∇d| p-4 ∇d • [∇d : ∇ 2 d]
with the property that f • d = 0.

Proof. For this, we consider d a 2-nematic solution, and compute its p-Laplacian:

∆ p d = ((d k,l d k,l ) p-2 2 d i,j ) ,j = (p -2)(d k,l d k,l ) p-4 2 d k,lj d k,l d i,j + |∇d| p-2 d i,jj = (p -2) |∇d| p-4 ∇d • [∇d : ∇ 2 d] + |∇d| p-2 ∆d, = (p -2)|∇d| p-4 ∇d • [∇d : ∇ 2 d] + |∇d| p-2 (-|∇d| 2 + f • d)d -|∇d| p-2 f = -(|∇d| p + |∇d| p-2 f • d)d -f , (57)
where f was defined in (56). Also we have that

(58) f • d = |∇d| p-2 f • d.
Indeed, the second term in (56) vanishes when multiplying by d as can be seen by differentiating the unitary condition d • d = 1. This yields (∇d) t d = 0, that is d i,j d i = 0, and then also

d k,lj d k,l d i,j d i = 0. In particular if f • d = 0 then also f • d = 0. Notice that if f • d = 0 then f • d = 0 so in this case λ p = 1 p |∇d| p .
Let us see, in the case of the 2-nematic two dimensional disk solution of Theorem 2, which is the p-nematic equation that this solution satisfies. We have the following result. and f ∈ L q (D) for all q < 4 p .

Proof. In this case, the solution is homogeneous (f = 0), so by Theorem 20 we have that where the integral is finite if q < 4/p, this complete the proof.

Summary and discussion

In this work we have obtained analytical solutions of the Frank-Oseen equations that are homogeneous or exact on a disk, and non-homogenoeus or approximate in a cylinder and a ball. The solutions possess infinite energy in the natural H 1 energy space, but they have instead finite energy in W 1,p spaces for some p < 2. We have shown that the exponent p is strictly smaller than 2 for the disk and cylinder cases, whereas for the sphere solution there are cases where the energy is finite only for p < 3/2. We have additionally defined a modification of the Frank-Oseen equations for allow infinite energy solutions as p-nematic liquid crystals. There is an interesting algebraic connection between 2-nematic and pnematic liquid crystals, whose regularity properties will be addressed in future studies. One of the main difficulties lies in the non-homogeneous loading term that arises on such a formulation, which we have shown it depends on the second order derivatives of the original solution. The solutions exhibit qualitative behavior with self-similar orbits or spirals in the horizontal plane and vortices in the vertical axis. This type of singularities naturally appear in nematic fluids. Finally, notice that some of the ansatz of the solutions were inspired in hydrodynamics orbits and singularities observed in proto-planet formation models. The rigorous connection between nematic fluids and these type of hydrodynamic or magneto-hydrodynamic models has not been yet established, and it will be matter of further research.

Figure 1 .

 1 Figure 1. Vector field and trajectories of the 2D disk solutions of Theorem 2 for n = 1 with R = 1 and c 1 = 3.84π.

Figure 2 .

 2 Figure 2. Trajectories of the 2D disk solutions of Theorem 2 for n = 1, n = 2 and n = 5 with R = 1 and c 1 = 3.84π.

Remark 4 .

 4 The family of solutions[START_REF] De Gennes | The physics of liquid crystals[END_REF] present a singularity of type log(r) at the origin. This means that the solution is in L 2 (D) but not in H 1 (D). We note that this singularity can be removed simply by considering the equation in a ring D\D with D = { ≤ |r| ≤ R}.

-

  sin a cos b cos a cos b -sin b   , e ⊥ (a, b) =   sin a sin b -cos a sin b -cos b   .

-

  sin a cn b cos a cn b -sn b   .

Figure 3 .

 3 Figure 3. Trajectories of the 3D cylindrical solutions for Theorems 5, 6 and 7 for R = 1 and H = tan(π/10). The coefficients used were n = 1, c 1 = 6π, c 2 = π for the first figure; n = 1, c 1 = 6π, c 2 = π/2 -arctan(H/R) for the second one; n = 1, c 1 = 3π, A = c 2 1 + n 2 , k = 0.87, c 2 = A/k for the third one.

Figure 4 .

 4 Figure 4. Protoplanetary disk HH-30 with jet, Taurus. Credit C. Burrows (STScI & ESA), the WFPC 2 Investigation Definition Team, and NASA [11].

Figure 5 .

 5 Figure 5. First component of the unitary field d from the cylinder case in the mid xy plane (left), sum in the z axis (center) and first component in the mid xz plane (right). The same three cases of Figure 3 are plotted in each row correspondingly.

Figure 6 .

 6 Figure 6. Streamlines for the spherical solutions of Theorems 11 and 13 in the unit sphere given by (42) and (49) respectively. Left: n = 1, c 1 = 3 and c 2 = -π/2. Right: n = 1, c 1 = 3, k = 0.87.

Figure 7 .

 7 Figure 7. First component of the unitary field d of the same spherical cases shown in Figure 6, first component of the solution (left), sum over z (center) and mid xz plane (right).

Figure 8 .

 8 Figure 8. Streamlines for spherical solutions in the spherical case of Theorem 13 for varying values of the parameter n. From left to right n = 1, 2, 5.

Theorem 21 .

 21 The 2-nematic solutions d ∈ W 1,p (D) for 1 ≤ p < 2 of Theorem 2 defined in the disk D = {(r, θ), 0 ≤ r ≤ R, 0 ≤ θ ≤ 2π} satisfy the non-homogeneous p -nematic equations of Definition[START_REF] Frank | On the theory of liquid crystals[END_REF] with right hand side given by(59) f (r, θ) = c 1 (p -2)2 p/2-2 (c 2 1 + n 2 ) p/2-1 r p/2 I ⊥ d(r, θ), (r, θ) ∈ D with f • d = 0, f = O 1 r p/2

2 I 2 |∇d| 2 = ∇λ 2 , 2 ,c 1

 222221 f = -(p -2)|∇d| p-4 ∇d • [∇d : ∇ 2 d].Changing variables we obtain∇ x,y d = 1 r ∇ r,θ d r cos θ r sin θ -sin θ cos θ = 1 ⊥ c 1 r d -nd r cos θ r sin θ -sin θ cos θ ,and we notice the following simplification (in tensor notation):∇d : ∇ 2 d = d k,l d k,lj = 1 2 (d k,l d k,l ) j = ∇ 1which allows us to compute∇d : ∇ 2 d = ∇λ 2 (x, y)Using that |∇d| p-4 = 2 p/2-2 λ p/2-2 we can finally compute the desired term:f = -(p -2)2 p/2-2 λ (p -2)2 p/2-2 (c 2 1 + n 2 ) p/2-1 r p/2 I ⊥ d,which naturally yields f • d = 0 and the required order. The regularity L q -regularity of f is given by

Acknowledgements A. O. was funded by ANID-Fondecyt 1201311, 1231404, CMM FB210005 Basal-ANID, FONDAP/15110009, Millennium Science Initiative Programs NCN19-161, ICN2021-004 and DO ANID Technology Center DO210001. This work was finished during the stay of A. O. at the Université Grenoble Alpes, Institut Fourier, CNRS UMR 5582. N. B. was funded by CMM FB210005 Basal-ANID and by ANID postdoctoral project number 3230326.