Anton A Bauhofer

Sebastian Krödel

Dr Jan Rys

Prof. C Osama R Bilal

Andrei Constantinescu

Chiara Daraio
email: daraio@caltech.edu

Daraio

Prof A Constantinescu

Harnessing Photochemical Shrinkage in Direct Laser Writing for Shape Morphing of Polymer Sheets

Keywords:

direct laser writing, residual stresses

Recently, residual stresses induced by polymerization were also used to produce structures with shape morphing potential. [11,12] Applications of programmable materials have been suggested in microrobotics, [13,14] biomedical devices, [15,16] and biomimetic systems. [6,17] However, most of the reported approaches require complex fabrication processes, including composites or multiple materials, or are limited to specific geometries or the use of soft materials. [18] Besides, most realizations of programmable materials have to date been demonstrated for macroscopic structures, in the centimeter and millimeter scales.

Here, we show a method to create microscale, programmable geometries by inducing residual stresses during direct laser writing (DLW) of photocurable polymers. DLW is a commonly used technique for the fabrication of microscale 3D structures. The technique is based on two-photon polymerization of photosensitive, liquid oligomers, generated by local exposure to a pulsed laser beam. [19][20][21] A byproduct of DLW is the presence of residual stresses, caused by incompatible strains due to nonuniform laser exposure, which result in undesirable shrinkage and distortions of the fabricated samples. The occurring deformations have been characterized in prior studies [22] and corrected through chemical alternations of photoresists [23] and geometrical precompensation of the expected deformations. [24] Here, we exploit these residual stresses to program thin, 2D sheets that reconfigure into 3D shapes. The reconfiguration is triggered by capillary forces, acting on the sheets during the drying process, after the fabrication steps of polymerization and developing. [25,26] The bending directions and angles in the resulting deformations are controlled by locally varying geometrical and writing parameters, opening a broad space for design. This approach allows for direct programmability of 3D-structured single-material sheets in a simple fabrication process, thus minimizing time and complexity in manufacturing.

The fabrication process and the bending mechanism are summarized in Figure 1. During DLW, a femtosecond laser is redirected and focused inside a droplet of photoresist, to induce local polymerization (Figure 1a). Due to the high laser intensity in the focal point, photoinitiator molecules in the photoresist are radicalized and trigger a chain reaction of the surrounding oligomers. In this process, oligomers are transformed into a dense polymer-network via photoinduced crosslinking Structures that change their shape in response to external stimuli unfold possibilities for more efficient and versatile production of 3D objects. Direct laser writing (DLW) is a technique based on two-photon polymerization that allows the fabrication of microstructures with complex 3D geometries. Here, it is shown that polymerization shrinkage in DLW can be utilized to create structures with locally controllable residual stresses that enable programmable, selfbending behavior. To demonstrate this concept, planar and 3D-structured sheets are preprogrammed to evolve into bio-inspired shapes (lotus flowers and shark skins). The fundamental mechanisms that control the self-bending behavior are identified and tested with microscale experiments. Based on the findings, an analytical model is introduced to quantitatively predict bending curvatures of the fabricated sheets. The proposed method enables simple fabrication of objects with complex geometries and precisely controllable shape morphing potential, while drastically reducing the required fabrication times for producing 3D, hierarchical microstructures over large areas in the order of square centimeters.

Programming shapes in planar materials enables the transformation of 2D sheets into controllable, 3D geometries. [1,2] This programmability expands the potential of available microand nanofabrication techniques and increases the fabrication throughput for objects with complex geometries [3] or origami structures. [4] Various physical effects have been exploited to produce structures that can change shape in a controlled manner, allowing for efficient fabrication of versatile geometries. [1] Among others, residual-stress-driven approaches have been explored to acquire self-bending behavior through gradients in thermal expansion, [5] anisotropic swelling, [6] differential lattice constants, [7,8] or buckling of prestrained substrates. [9,10] (Figure 1b). The formation of covalent bonding between the oligomers locally reduces the average atomic distance in the resist, [27] leading to an increased density at the molecular level. At the macroscopic level, this induces shrinkage and distortions of the polymerized volumes (Figure 1c). By scanning the laser focal point along a prescribed trajectory, regions that experience shrinkage can be formed with almost arbitrarily complex geometries. These regions can then be utilized to evoke shapechanging behavior. Laser writing a sheet of parallel adjacent voxel lines leads to crosslinking of the polymerized regions, causing collective in-plane shrinkage (Figure 1d). During consecutive writing of multiple stacked layers, interconnectivity of the layers inhibits further in-plane shrinkage and results in the emergence of residual stresses that create an out-of-plane bending moment (Figure 1e). After laser writing, a developing process dissolves and washes off the remaining resist. During drying, the surface tension of the washing liquid outweighs adhesive forces between the structures and the underlying substrates, triggering delamination of the sheets from the substrates and out-of-plane bending. We fabricated and analyzed sheets of stacked layers with a thickness of t = 0.5 µm and lateral dimensions of 250 µm × 250 µm (Figure 1f). The resulting bending curvatures of the sheets are strongly dependent on the laser writing trajectories and the number of stacked layers, which dictate the magnitude of the residual stresses as well as the bending rigidity of the sheets. Stress concentrations within the sheets cause delamination to begin from the sheet edges and then propagate inward (see Supporting Information).

To understand the fundamental mechanisms of self-bending induced by photopolymerization, we designed controlled experiments to identify the driving physical effects. We determined that i) capillary action, ii) residual stresses, and iii) adhesive forces are the main influence factors controlling the resulting bending curvatures and geometries (Figure 2a). First, we investigated the influence of capillary action on the bending, which occurs after developing of the polymer, by applying different drying procedures and rinsing liquids to separate samples. To test the influence of different surface energies of the rinsing liquid, we immersed developed samples in water and isopropanol (IPA) and measured the receding contact angle of the liquids on the sheets (Figure 2b). When the samples were dried in water, we observed evaporation with a pinned contact line at the sheet edges and a receding contact angle of θ R,water = 40.5° ± 0.1°. Due to the high surface energy of the liquid, the tested samples were completely detached from the substrates. During drying of samples in IPA, we observed no pinned contact line, due to the low surface energy of IPA, and a constant contact angle of θ R,IPA = 16° ± 1°. By using IPA as drying liquid, we achieved controlled repeatable detachment and self-bending. We also performed control experiments using critical point drying (CPD), where no capillary actions occur. In this case, the samples did not show any significant bending, demonstrating the relevant role played by capillary action on the programmability of the final shapes.

Next, we analyzed the residual stresses that form during photopolymerization. For the chosen DLW parameters (see Experimental Section), we identified voxel lines with ellipsoidal cross-sections with a width and height of 0.37 ± 0.03 µm and 2.07 ± 0.20 µm and a longitudinal strain of ε l = 16.3 ± 1.1% that occurs during solidification of the photoresist (see Supporting Information). Crosslinking between consecutively written voxel lines restrains shortening of the lines, resulting in a reduced in-plane macroscopic shrinkage within sheets. We used optical microscopy to observe photochemical shrinkage in-situ during DLW and evaluated the macroscopic shrinkage of sheets with scanning electron microscopy (SEM). For ten evaluated, unrestricted single-layer sheets with parallel voxel lines, a mean polymerization shrinkage-induced strain of ε p = 3.7 ± 0.5% was found. This value agrees with prior literature. [28] The consequent residual stresses that occur during DLW are calculated as σ r = (1 -ν) -1 × ε p × E, where E and ν denote the Young's modulus and the Poisson's ratio of the polymerized material. To characterize the mechanical properties of the polymer sheets, we applied nanoindentation (see Supporting Information). We obtained a Young's Modulus of E = 1.86 ± 0.15 GPa in accordance with prior literature, [29] while assuming a Poisson's ratio of 0.49 υ =

. [29] For measuring the lateral adhesive forces between the polymer sheets and adjacent glass substrates, we used a micromechanical testing system (see Experimental Section). We pushed a laser written cube on the glass surface and measured the resulting force with a capacitive force sensor (Figure 2c We develop a mechanical model to inform the experimental observations and quantitatively predict delamination and self-bending angle of multilayered sheets for design. In the model, we account for the fact that detachment and bending angle are both highly sensitive to the number of stacked layers n (Figure 1f). From experimental observations, we consider two distinct regimes: i) sheets with a low number of layers (n ≤ n D ) that undergo no delamination, and ii) sheets with a high number of layers (n > n D ) that experience delamination and subsequent self-bending into specific curvatures. n D is thereby introduced as a critical number to locate the point of delamination.

To determine the critical number n D , we derive a delamination condition that balances the in-plane forces acting in tangent direction at the delaminating edge (Figure 2a). Line loads in an element with infinitesimal width dx are calculated to capture the experimentally measured effects. The line load imposed by capillary action during drying is computed as q s = cos θ R × γ lv , where we assume the surface tension between the liquid and the vapor phase of IPA to be γ lv, IPA = 21.0 × 10 -3 N m -1 . The impact of residual stresses is computed by integration of the residual stresses over the total sheet height d

rr ∫ σ = × qz t nt
, where the thickness of the stacked layers is denoted as t. We assume a uniform stress distribution in the cross-sections of the stacked layers (ε i = ε p , i > 1), while the first layer (i = 1) is presumed to experience unconstrained shrinkage. The total compressive line force at the delaminating edge amounts to the sum of the line loads due to the two effects q c = q s + q r . The maximum adhesion between the line element and the substrate is calculated as q a,max = τ max × w, where w denotes the lateral sheet dimension. Delamination is predicted if the contractile forces outweigh the adhesion of the line element q c > q a . Based on the experimental results, the critical number is computed as n D = 3.30 ± 0.86 (Figure 2d). The self-bending curvatures of the sheets after delamination are predicted with a simplified model based on multilayer beam bending theory. Sheets in a cylindrical bending configuration are considered in a state of minimal stresses. Curvatures normal to the cylindrical bending are assumed to be zero, thus reducing the model to a 1D system that can be approximated with beam theory. All stresses in an element of the bent beam portion are represented by n, tensile forces P i , and bending moments M i for the distinct layers (see Supporting Information). [30,31] In the final sheet configuration, no external forces or moments act on the bent region of the sheets, leading to force and moment equilibria that only contain the introduced loads P i and M i . The constitutive equations for the layers are

M EI i ii ρ =
, where E i and I i represent the Young's Moduli and the moments of inertia, while ρ is referred to as the radius of curvature and is defined as the reciprocal value of the curvature 1 ρ κ = . In addition, we introduce a compatibility condition for the strains in adjacent layers We assume that cross-sections which are plane and perpendicular to the neutral axis remain in that condition and that out-of-plane tensions related to edge effects are negligible. Furthermore, the sheet thickness is assumed to be much smaller than the lateral sheet dimensions

 tw i n ∑   
 , while the distinct layers are approximated to be homogeneous and isotropic. The resulting closed-form solution (see Supporting Information) is

6 > 1 1 3 ∑ ∑ ∑ κ εε ( ) ( ) = - = = tt t i ji n i n ji j i i n (2)
To achieve controlled adhesion, the first written layer is partly immersed in the substrate. Therefore, we assume half of the thickness for the first layer

1 2 1 =       tt (see Supporting Infor- mation)
. For all proceeding layers, the thickness is assumed to be constant (t i = t, i > 1). With these assumptions and the introduced delamination condition, the solution expressed in Equation ( 2) can be simplified to 0, 24

1 81 26 1 ,> D p 32 D nn t n nn n nn κ ε = ≤ × - -+ -      (3)
We compare the predicted curvatures from the model to the experimentally determined curvatures of laser written sheets (Figure 2f). The model and the experimental results agree for both the nondelaminated and the self-bent regions. We determined the bending angles β of the bent sheets with SEM imaging, by measuring the tangential slopes of the free edges (Figure 2e). The viewing angle δ was corrected with cos β β δ = ′ , where β′ denotes the projected bending angle in the SEM image. We computed the curvatures from the bending angles as

κ β = l
, where the arc length is denoted as l. The presented experimental approach, paired with our predictive analytical tool, enables the direct fabrication of microscale shape morphing structures with designed bending curvatures and shapes. This method can be used to produce 3D structures at a fraction of the time and fabrication costs compared to their realization by DLW (without shape morphing). To demonstrate the potential of our approach, we fabricated self-bending sheets that evolved into lotus-like microcontainers upon drying (Figure 3a). We controlled the bending curvatures by varying the number of stacked layers. Positive Gaussian curvatures were achieved by stacking of layers with differently oriented voxel lines. While the microcontainers could also be fabricated by conventional DLW, laser writing time was reduced by a factor of ≈4 when fabricated using self-bending sheets (the sheets required a writing time of T S = 5.5 ± 0.5 s, whereas the direct writing of the final, spherical geometry took T o = 19.0 ± 0.5 s, see Supporting Information). A potential application for the self-bending geometry can be found in drug delivery, where the approach could be applied to biodegradable materials that exhibit photochemical shrinkage. [16,32,33] Our method can also be used to create programmable anisotropic and inhomogeneous surfaces that replicate functionalities found in nature (Figure 3b). For example, shark skin topology is known to enhance swimming speed and agility of the fish, by creating a leading-edge vortex. [34] We fabricated a biomimetic shark skin in the form of a 3D-structured sheet with structured microfins. To locally control bending, we increased the sheet thickness in selected locations (see inset of Figure 3c). These variations in thickness enable localized selfbending and lead to the formation of structures that mimic the geometry and feature sizes of natural shark skins. Using DLW, surfaces in the cm-range can be produced in a single production cycle (see Supporting Information).

In conclusion, we show that DLW can be used to create microstructures with programmable self-bending behavior. All driving physical effects for the acquired shape transformations have been identified and evaluated using ad-hoc microscale experiments. We present a mechanical model that predicts bending curvatures based on the measured parameters. Due to the flexibility of the fabrication process, our approach opens prospects for a fast and simple creation of shape morphing structures with complex geometries.

Experimental Section

Photolithography: The self-bending microstructures were fabricated with the commercial DLW system Photonic Professional GT by Nanoscribe GmbH. The photoresist IP-Dip (Nanoscribe GmbH, Eggenstein-Leopoldshafen, Germany) and a 25× objective with a numerical aperture of 0.8 and a working distance of 390 µm were used in galvo scan mode in dip-in laser lithography configuration. Laser power (LP = 16 mW), scan speed (SS = 55 000 µm s -1 ), and galvo settling time (GT = 2 ms) were kept constant during writing, while slicing and hatching distances were set to 0.5 and 0.3 µm. The implemented parameters and materials serve as a benchmark whereby qualitatively equivalent results would be expected for the use of other resists. The observed effects could be potentially maximized with resists that feature stronger polymerization shrinkage. [35] Multilayered sheets of parallel voxel lines were made with the dimensions 250 µm × 250 µm × (n -0.5) × 0.5 µm with n being the number of stacked layers. The laser written structures were developed with the solvent mr-Dev 600 (Micro Resist Technologies GmbH, Berlin, Germany) for 12 min, rinsed in IPA for 3 min, and air-dried in horizontal orientation in steady air. The receding contact angle of IPA on the samples during evaporation was measured with a DSA100 drop shape analyzer (Krüss Ltd., Hamburg, Germany). To verify that capillary forces act as trigger for the self-bending behavior of the laser written sheets, CPD, during which no capillary actions occur, was performed in comparison to air drying. No self-bending curvatures were observed in sheets that were subjected to CPD.

Micromechanical Testing: We applied a micromechanical compressions sense system (FemtoTools AG, Zurich, Switzerland) for measuring the adhesive forces between laser written structures and uncoated glass substrates. A FT-L1000 sensing probe with a force range of ± 1000 µN and a resolution of 0.05 µN was used with a 50 µm × 50 µm tip to measure the occurring forces during lateral movement of blocks of polymerized photoresist with the dimensions 25 µm × 25 µm × 25 µm. The material properties of the constituent polymer film were measured using a nanoindenter (TriboIndenter, Hysitron Inc., USA) equipped with a 5 µm, 60° conical tip. All measurements were carried out in a displacement-controlled mode with an indentation depth of 2.0 µm. Each measurement consisted of a loading function including a loading (10 s) and an unloading (10 s) segment. The material properties of five individual samples were determined on four different positions around the center of the film. The reduced Young's modulus was determined by the Oliver and Pharr method. [36]

Figure 1 .

 1 Figure 1. Direct laser writing (DLW) of prestrained sheets prone to self-bending. a) An 80 MHz femtosecond laser is redirected and focused with galvo mirrors and an objective into a droplet of photoresist to cause local photopolymerization in the focal point. b) Representative radical polymerization reaction of oligomers with acrylate end-groups. c) The chemical reactions at the molecular level induce a macroscopic volumetric shrinkage, denoted as polymerization shrinkage. d) Layers of polymerized material are created by laser writing of parallel voxel lines. Crosslinking of adjacent voxel lines causes shrinkage in all lateral directions. e) Successive shrinkage of stacked layers induces a bending moment that enables out-of-plane deflections. f) SEM images of different self-bent sheets containing 1 to 10 (from left to right) stacked layers of parallel voxel lines that were successively written from bottom to top. Scale bar: 100 µm.

  ). A maximum lateral adhesion tension of cos 623was determined, where F m,max = 321.7 ± 61.9 µN denotes the maximum measured force, α = 15° the mounting angle of the sensor, and A m = 625 µN 2 the surface area of the test cube.

Figure 2 .

 2 Figure 2. Self-bending mechanism of laser-written sheets, detaching from the substrate. The shaded areas in panels c-f indicate the standard deviations computed from the measurement results with Gaussian error propagation. a) Schematic diagram showing the lateral forces acting along the border of detachment on a line element with infinitesimal width dx. b) Measured receding contact angles of water and IPA on a substrate produced by DLW. The inset shows the contact angle of a water droplet. c) Measurements of the lateral adhesive forces between the photopolymer and the glass substrate. The inset depicts the experimental setup. d) The experimentally determined contractile line loads from the capillary action and the residual stresses are compared to the maximum adhesive line force q a,max = τ max × w. Delamination is expected for sheets with a number of layers above the critical number n D = 3.30, when the contractile forces overcome the maximum adhesion. e) Self-bent sheet with an angle of curvature β, a radius of curvature ρ, a width w, and an arc length l. Scale bar: 100 µm. f) Comparison of the analytical predictions with experimental data for the final bending curvature obtained with samples composed of different numbers of layers. The error bars indicate the standard deviation determined from measurements on 10 samples per data point.

Acknowledgements

This work was partially funded by the Swiss National Science Foundation through grant "MechNanoTruss-Mechanical response of polymer nanotruss scaffolds" (No. 164375).

Conflict of Interest

The authors declare no conflict of interest.


Krister M Swenson
email: swenson@lirmm.fr

Afif Elghraoui
email: aelghraoui@sdsu.edu

Faramarz Valafar
email: faramarz@sdsu.edu

Siavash Mirarab
email: smirarab@ucsd.edu

Mathias Weller
email: mathias.weller@cnrs.fr

Elghraoui, Valafar, Mirarab Weller Swenson

Quantifying Hierarchical Conflicts in Homology Statements

Keywords: Homology, Syntenic Block, T-Star Packing, Assignment Problem

A fundamental step in any comparative whole genome analysis is the annotation of homology relationships between segments of the genomes. Traditionally, this annotation has been based on coding segments, where orthologous genes are inferred and then syntenic blocks are computed by agglomerating sets of homologous genes into homologous regions. More recently, whole genomes, including intergenic regions, are being aligned de novo as whole genome alignments (WGA). In this article we develop a test to measure to what extent sets of homology relationships given by two different software are hierarchically related to one another, where matched segments from one software may contain matched segments from the other and vice versa. Such a test should be used as a sanity check for an agglomerative syntenic block software, and provides a mapping between the blocks that can be used for further downstream analyses. We show that, in practice, it is rare that two collections of homology relationships are perfectly hierarchically related. Therefore we present an optimization problem to measure how far they are from being so. We show that this problem, which is a generalization of the assignment problem, is NP-Hard and give a heuristic solution and implementation. We apply our distance measure to data from the Alignathon competition, as well as to Mycobacterium tuberculosis, showing that many factors affect how hierarchically related two collections are, including sensitivities to guide trees and the use or omission of an outgroup. These findings inform practitioners on the pitfalls of homology relationship inference, and can inform development of more robust inference tools.

Introduction

The increasing ease of whole genome sequencing and assembly has opened a new era of comparative genomics. With the data available today, not only can the phylogenetic histories of all the genes between a set of genomes be analyzed, but also the interaction between these genes, linking gene regulation and function to the positions of groups of genes. These analyses require a reliable grouping of homologous genomic segments from the multiple genomes in question.

Thus, the inference of sets of homologous genomic segments is of fundamental importance. Such a homology statement comes in the form of a set of genomic segments that contains at least one, but potentially multiple, segments from several genomes. Each pair of segments from the set shares common ancestry over some proportion of their intervals, which varies depending on the scale and level of precision required by the application.

The most basic segment on which statements are made has traditionally been the gene, detected through either manual or automatic means. The number of tools designed to infer homology relationships between annotated genes has grown, provoking the formation of the Quest for Orthologs (QfO) consortium dedicated to the evaluation and comparison of these tools [START_REF] Glover | Advances and applications in the quest for orthologs[END_REF].

General genomic intervals, that can contain both coding and noncoding positions, have also been used as homology statements. In this case, researchers have considered bidirectional best hits as evidence for orthology [START_REF] Mural | A comparison of wholegenome shotgun-derived mouse chromosome 16 and the human genome[END_REF][START_REF] Waterston | Initial sequencing and comparative analysis of the mouse genome[END_REF]. More recently, "whole genome alignment" methods partition entire genomes into blocks that can be aligned into multiple sequence alignments (MSAs), de novo, with no special input from the user. The Alignathon collaborative competition was developed to evaluate and compare these methods [START_REF] Earl | Alignathon: a competitive assessment of whole-genome alignment methods[END_REF].

Study of the large scale changes between genomes has inspired a more vague notion of homology between genome segments. Even before the discovery of the double helix, groups were studying homology of large segments of genomes from the salivary glands of drosophila [START_REF] Sturtevant | The homologies of the chromosome elements in the genus drosophila[END_REF]. More precise lengths of roughly conserved chromosomal segments began to be studied using linkage maps [START_REF] Nadeau | Lengths of chromosomal segments conserved since divergence of man and mouse[END_REF]. In the postgenomic era, basic homology segments are agglomerated into syntenic blocks, possibly separated by micro-rearrangements. GRIMM-synteny was developed for the study of large scale chromosomal changes, in response to the whole genome sequencing efforts in human and mouse [START_REF] Pevzner | Genome rearrangements in mammalian evolution: lessons from human and mouse genomes[END_REF]. Since then, many syntenic block inference tools have been introduced but, despite twenty years of development, a unified definition of syntenic block has yet to be found. Indeed, most tools rely on operational definitions rather than biological or mathematical ones [START_REF] Ghiurcuta | Evaluating synteny for improved comparative studies[END_REF][START_REF] Svetlitsky | Discovery of multi-operon colinear syntenic blocks in microbial genomes[END_REF].

There are several tests used for comparing and evaluating homology statements. For orthology statements between coding sequences, the QfO project has established tests that: 1) compare trees inferred from orthologous families to agreed-upon species trees, 2) compare the subsets of orthologs from curated gene families, and 3) use consistency in gene ontology annotation [START_REF] Altenhoff | Standardized benchmarking in the quest for orthologs[END_REF]. For statements between syntenic blocks, comparisons of inclusion and exclusion of segments between methods have been done, and blocks from a single method have been compared to known gene clusters [START_REF] Marcet-Houben | EvolClust: automated inference of evolutionary conserved gene clusters in eukaryotes[END_REF]. Other sources of ground truth, such as RNA-seq data, have been used to confirm co-regulation between genes occurring in a proposed block [START_REF] Winter | Finding approximate gene clusters with gecko 3[END_REF].

For general statements (on both coding and noncoding DNA), the Alignathon competition used three different measures [START_REF] Earl | Alignathon: a competitive assessment of whole-genome alignment methods[END_REF]. If homology statements are given as a set of (potentially gapped) equal-length segments from several genomes, then each homologous pair of positions between two genomes as given by one method, can be queried in another method. The number of such shared positions is a measure of similarity and, when one of the methods is taken as ground truth, the number of shared positions can be used to measure precision and recall. The mafComparator tool estimates these values by sampling positions [START_REF] Earl | Alignathon: a competitive assessment of whole-genome alignment methods[END_REF]. For sets of aligned homology statements (i.e. MSAs), probabilistic sampling-based alignment reliability (PSAR) was used to assess each aligned column [START_REF] Kim | PSAR: measuring multiple sequence alignment reliability by probabilistic sampling[END_REF]. PSAR fixes all rows of the alignment but one, and samples from the many ways to align that row within the fixed alignment. After this is repeated for each row, an alignment reliability score for each pair of positions in a column can be assigned. When aligned homology statements are augmented with a phylogeny, another statistical test called StatSigMA can be used [START_REF] Prakash | Measuring the accuracy of genome-size multiple alignments[END_REF]. For each edge of the phylogeny, the rows of the alignment are split into two alignments. The two alignments are then tested for exhibiting "unrelated behaviour" using Karlin-Altschul log likelihood scores. If the test for all branches passes, then the homology statement is validated.

For homology statements that come in the form of syntenic blocks, Ghiurcuta and Moret outline some necessary conditions for a valid agglomeration of homologous units into such blocks [START_REF] Ghiurcuta | Evaluating synteny for improved comparative studies[END_REF].

There exists very few methods that compare homology statements in the form of sets of genomic segments, unmarried to connotations of orthologous genes, and independent of multiple sequence alignments. To our knowledge, the Jaccard distance (e.g. as computed by mafComparator) applied to pairwise homology statements, is the only known comparison that falls into this category.

In this article we introduce a simple definition of homology block (Section 2.1) and formally characterize the conditions under which a set of homology blocks are valid (Section 2.2). We show what it means for collections of blocks to be hierarchically related and use this to develop a method for measuring disagreement between two different collections of blocks. In Section 2.3 we show a necessary condition for two collections of blocks to be in a hierarchical relationship (in the form of Lemma 1), based on a graph representing the overlap between the sets. For different parts of the genomes in question, our test allows for the collections of blocks to be hierarchically related in both ways; in some parts of the genomes the first set could be more general than the second, while in other parts of the genomes the opposite can be true. We introduce an optimization problem, called Minimum Deletion into Disjoint Stars (MDDS), which gives a lower bound on the number of positions that must be ignored so that the two collections of blocks could be related through a hierarchical relationship. Not only does a solution to MDDS give a measure as to the degree of hierarchical dissonance between two collections, but it serves as an unambiguous mapping between the blocks of the two. This mapping could be used for further downstream comparisons, an illustration of which is shown in Appendix C.

We show that the MDDS problem is NP-Hard, before presenting a polynomial time heuristic based on an exact algorithm for solving MDDS on a tree. In Section 4 we define the homology discordance ratio and use this measure as a distance between block collections built on Alignathon data and on a set of 94 Mycobacterium tuberculosis isolates. On the tuberculosis strains we study the relationship between blocks built using an outgroup or no outgroup, using annotations or no annotations, using maf2synteny to agglomerate or not, as well as study the effect of the guide tree on Cactus MSA blocks. On the simulated data from the Alignathon project, we highlight differences between our method and the Jaccard distance (as computed with mafComparator).

For the entirety of the article we focus on the general case of homology statements, although most of the discussion also applies to the restricted case of orthology statements.

Methodological Foundations

Overlapping Homology Statements and the Block Graph

We use g i [k..ℓ] to denote a segment between positions k and ℓ of genome g i and we let T denote the universe of all segments over all possible genomes and position-pairs. Define the overlap op(s 1 , s 2 ) between two segments s 1 and s 2 as the number of positions where they overlap in the same genome. For example op(g 1 [1..5], g 1 [3..9]) = 3 but op(g 1 [1..5], g 2 [3..9]) = 0.

Definition 1 (homology statement block). A homology statement block (called a block for short) B is a set of segments B ⊂ T such that all pairs of segments in B have zero overlap.

The right panel of Fig. 1 Before discussing the semantic interpretation of a homology statement block, we first introduce a graph that represents the overlap between blocks. The overlap op(B1, B2) = s1,s2∈B1×B2 op(s 1 , s 2 ) between blocks B1 and B2 is the total overlap between all pairs of segments in the two. A collection of blocks B = {B1, B2, . . .} is considered to be clean if the overlap between any pairs of blocks in B is zero. Both collections depicted in Fig. 1 are clean.

For two collections of blocks A and B, we build a bipartite block graph BG(A, B) where there is an edge between A and B for any A ∈ A and B ∈ B if Segments are depicted with brackets and lined up according to their positions on the chromosomes. They are labeled by their tuple (when space permits) and the block to which they belong. The configuration of positive and negative witness pairs shows that B generalizes A. Some of the genome positions are highlighted with boxes, and two such positions are connected by a solid line if they appear as a positive witness in B, and that line is bold if they are also in A. The dashed line represents one of the (many) negative homology witnesses between A1 and A2 that are negative witness pairs for A but not for B. To the left, the graph BG(A, B) appears with edges labeled by overlap length in gray. All of the connected components are stars. and only if blocks A and B overlap (i.e. op(A, B) > 0). Thus, the block graph BG(A, A) for a clean collection A is a perfect matching. E(G) is the set of edges of the graph G. We associate to each edge AB ∈ E(G) a weight function ω : E(G) → N such that ω(AB) = op(A, B). The left side of Fig. 1 shows the block graph for the collections to the right.

Homology Witnesses and Block Hierarchies

A homology block can be interpreted as a positive and negative statement of homology (i.e. statements about common ancestry). On the positive side, the block {g 1 [1..5], g 2 [11..16]} says that positions 1 through 5 in genome g 1 are somehow homologous to positions 11 through 16 in genome g 2 (in this case the segments are not the same length, so we assume that each position from the segment g 1 [1..5] is homologous to either a position in g 2 [11..16] or to none other in g 2 ).

On the negative side, the block could be interpreted as saying that no position in g 1 [1..5] is homologous to any other position in g 1 or any other position in g 2 , outside of g 2 [11..16].

In this section we suppose that we know the truth about the ancestral relationships between the base-pair positions of the genomes in question. With this supposed knowledge, we can categorize pairs of homology witness positions as positive or negative, depending on their evolutionary relationship. Using these relationships, we define properties that a valid collection of homology blocks must respect. These definitions are extended to encompass hierarchical relationships between collections of blocks.

Consider any pair of positions g i [x] and g j [y] from genomes g i and g j . This pair is called a positive homology witness if the two positions descend from the same ancestral position, otherwise the pair is called a negative homology witness (positive homology witnesses represent pairs of positions that are typically called "homologous" positions). Note that the true relationship between positions is unknown, yet it imposes constraints on what we consider a valid collection of blocks according to the following definition.

Consider any position-pair

(g i [x], g j [y]) such that g i [x] is contained in a seg- ment from a block B in a collection B. If (g i [x], g j [y]
) is a positive homology witness, then either 1. g j [y] appears in B, and we say that the pair is a positive witness in B, or 2. g j [y] appears in no block of B.

By this definition no position-pair (g i [x], g j [y]) with g i [x] and g j [y] in different blocks of B, can be a positive homology witness and, since all position-pairs are either positive or negative homology witnesses, (g i [x], g j [y]) must be a negative homology witness. Any position-pair (g i [x], g j [y]), where g i [x] and g j [y] are in different blocks of B or in no block of B, is called a negative witness for B.

Note that, for a clean collection B, no position-pair can be both a positive and negative witness in/for B. There may also be position-pairs that are neither positive nor negative witnesses in/for B, such as those pairs that have one position contained in a block of B and the other outside all blocks of B. Finally, note that not all position-pairs appearing between segments in a homology block need necessarily be positive homology witnesses.

Positive witness pairs limit what can exist in two different blocks; a block containing one position of a positive homology witness imposes the constraint that the other position must either be in the same block, or in no block. On the other hand, we will see in the following that negative homology witnesses existing between two different blocks in a collection enforce constraints on the hierarchical relationships that this collection can have with another block collection.

Consider the collections of blocks in Fig. 1 and note that whenever a positive homology witness is a positive witness in A, it must also be a positive witness in B, whereas not all positive witnesses in B exist in A. Conversely, every negative witness for B is also a negative witness for A. In this sense, the blocks of B are "more general" than the blocks of A. This motivates the following definition, for which we focus on subcollections of blocks A ′ ⊆ A and B ′ ⊆ B.

Definition 2 (generalization).

A clean (sub)collection of blocks B ′ generalizes a clean (sub)collection A ′ if and only if every positive witness in A ′ is also a positive witness in B ′ and every negative witness for B ′ is also a negative witness for A ′ .

Note that any clean collection generalizes itself.

While some subcollections of B may generalize subcollections of A, other subcollections of A may generalize subcollections of B.

Partition them A = A 1 ∪ A 2 ∪ • • • ∪ A k and B = B 1 ∪ B 2 ∪ • • • ∪ B k according to the connected components of BG(A, B) (e.g. A 1 ∪ B 1 is the set of vertices in the first connected component).
Definition 3 (hierarchical). We say that A and B have a hierarchical relationship if and only if

A i generalizes B i , or B i generalizes A i , for 1 ≤ i ≤ k.
The existence of hierarchical relationships between collections of blocks are interesting to us for at least two reasons. Consider two block inference methods, MethodA and MethodB, producing different clean collections of blocks A and B respectively. If MethodB is meant to agglomerate blocks from MethodA, then we would expect B to generalize A. This is useful for the verification of agglomeration methods, and as a sanity check for the practitioner. In this case, if MethodB also trims spurious blocks or segments from MethodA, B may not generalize A, but A and B would still be hierarchically related. Another reason for interest in the hierarchical relationship may be that, if B generalizes A, then we can define a mapping from each block A ∈ A to a block B ∈ B. This mapping can be used for further comparisons between the collections. The refinement of orthology assignments, as illustrated in Appendix C is an example of one such comparison.

Relating Block Hierarchy to Stars in the Block Graph

While simple hierarchical relationships are easy to detect, real-world data are not so well behaved, and require a formalism to measure the extent to which a relationship is hierarchical. The types of connected components in the block graph give us insight into collections that cannot have a hierarchical relationship. Lemma 1. Let A and B be clean collections of blocks such that B generalizes A and BG(A, B) is connected. Then, all vertices of A have degree one in BG(A, B), that is, BG(A, B) is a star with center in B.

Proof. Let A be a block in A, and assume that it has at least two distinct neighbors B1, B2 ∈ B in BG(A, B), that is, both B1 and B2 overlap A. Thus, there are positions g i [x] and g j [y] appearing in A such that g i [x] appears in B1 and g j [y] appears in B2. Since B is clean, we also know that these positions are distinct. Since (g i [x], g j [y]) appears in different blocks in B, we know that it is a negative homology witness for B. However, since g i [x] and g j [y] appear in the same block in A, the pair is not a negative witness for A. Thus, (g i [x], g j [y]) is a negative witness for B but not for A, contradicting the fact that B generalizes A.

⊓ ⊔

We say that a graph is hierarchical if it is a collection of vertex-disjoint stars, that is, if no component has two vertices of degree greater than one. It is easy to check if a graph meets this criterion. Note that the condition of a graph BG(A, B) being hierarchical is necessary for A and B to have a hierarchical relationship, but it is not sufficient. Also note that, Lemma 1 outlines a property on each individual connected component, allowing some parts of A to generalize parts

A1 A1 A1 A1 A1 A1
Fig. 2: Two subsets of blocks from collections A and B that are not hierarchically related. The top panel depicts a subset of blocks for part of the genome g 1 , but not the other genomes. Positions α and β form a negative homology witness for A, but not for B, while β and γ form a negative homology witness for B, but not for A. This contradicts properties of a hierarchy and, therefore, yields the nonstar topology to the left. The bottom panel depicts a subset of blocks for part of the genomes g 1 and g 2 . These segments contradict properties of a hierarchy in a different way, and yield the non-star topology to the left. This kind of scenario would arise when orthologs are matched in one way from MethodA, and in another way for MethodB. Note that even if the block A5 was not in collection A, the contradiction still holds and the connected component is not a star.

of B, while allowing other parts of B to generalize parts of A. Thus, a natural corollary to Lemma 1 is that if A and B are hierarchically related, then BG(A, B) is hierarchical.

For a graph that is not a collection of stars, one may want to measure to what degree it deviates from being so. Lemma 1 inspires the search for star packings on G = BG(A, B).

Input: a bipartite graph G with weight function ω : E(G) → N Output: E ′ ⊆ E(G) such that the subgraph of G formed by the edge set (E(B) \ E ′ ) is a collection of vertex-disjoint stars Measure: e∈E ′ ω(e)

Minimum Deletion into Disjoint Stars (MDDS)

A solution to MDDS gives a lower bound on the number of overlapping positions that must be ignored so that A and B can be hierarchically related. For example, Fig. 2 shows two connected components that are not stars. Consider the graph from the upper panel and assume that in the non-depicted genomes (i.e. g i for i > 1) there is no overlap of the segments of A1 with those of B2, or of segments of A2 with those of B1. Then, the solution to MDDS on this component would result from the removal of the edge between B1 and A2, since this edge has the minimum-size overlap. In Section 4.1 we highlight the differences between our MDDS method, and the Jaccard similarity index used in the Alignathon.

Algorithms

In this section we show that the Minimum Deletion into Disjoint Stars problem is NP-Hard, and then present a practical heuristic based on solving MDDS optimally on a tree. For simplicity and without loss of generality, we will assume that all block graphs are connected.

Other generalizations of the assignment problem, similar to MDDS, have been studied for decades, the closest of which has most recently been called the T -star packing problem [START_REF] Babenko | New exact and approximation algorithms for the star packing problem in undirected graphs[END_REF]. This problem asks for a star packing where the size of the star is limited by an input parameter T . When the measure is the number (or weight) of edges, Hell and Kirkpatrick show that the T -star packing is NP-Hard by reduction to the version of the problem that asks for a decomposition of a given graph into subgraphs isomorphic to the star with T edges [START_REF] Hell | Packings by complete bipartite graphs[END_REF]. Since the only difference between MDDS and the edge-weighted T -star packing problem is the parameter T , it is tempting to adapt the same series of reductions to MDDS by setting T to the maximum degree over all vertices in the graph. This approach is not clearly feasible, however, since the reduction from 3-Dimensional Matching to the decomposition version of the problem creates vertices of degree higher than T [START_REF] Kirkpatrick | On the complexity of general graph factor problems[END_REF].

Babenko and Gusakov give a 9 4 T T +1 approximation algorithm for the T -star packing problem based on a reduction to the max-network flow problem [START_REF] Babenko | New exact and approximation algorithms for the star packing problem in undirected graphs[END_REF]. We could use this elaborate approximation algorithm by fixing T to the maximum degree of the input graph, but we choose instead to implement the much simpler heuristic presented in Section 3.2.

NP-Hardness of MDDS

We will show that the decision version of MDDS is NP-hard by reducing the well-known 3-SAT problem to it. Our construction uses similar techniques as the NP-hardness proof of the Transitivity Edge Deletion problem [START_REF] Weller | On making directed graphs transitive[END_REF].

Construction 1 (see Fig. 3) Let φ be an instance of 3-SAT with variables X := {x 1 , x 2 , . . . , x n } and clauses C := {C 1 , C 2 , . . . , C m } such that each clause contains exactly three literals. For each variable x i , let n i denote the number of clauses that contain x i or ¬x i and let γ 0 i , γ 1 i , . . . , γ ni-1 i be any sequence of these clauses. We construct an edge-weighted graph (G, ω) as follows: 1. For each variable

x i create a cycle Q i containing 6n i vertices v 0 i , v 1 i , . . . , v 6ni-1 i
and give all edges weight m. , that is, it is the 4 th clause containing x 1 , the 2 nd clause containing x 3 and the 6 th clause containing x 8 . A truth-assignment setting x 1 to TRUE and x 3 and x 8 to FALSE corresponds to the star cover indicated by gray highlights. Note that taking the edge between v 18 1 and C instead of the edge between v 32 8 and C corresponds to satisfying the clause C by x 1 instead of ¬x 8 .

3.

For each i, j let ℓ be such that γ j i = C ℓ and, if x i occurs non-negated in C ℓ , then add the edge {v 6j i , u ℓ } with weight 1, otherwise add the edge {v 6j+2 i , u ℓ } with weight 1. Note that the image of ω is {1, m}, the total weight of all edges is 18m 2 + 3m, and G is bipartite, since any edge from part 3 of the construction, connecting a (variable) cycle to a (clause) vertex u ℓ , connects to an even numbered vertex in the cycle.

Besides NP-hardness, our reduction implies exponential lower bounds assuming widely believed complexity-theoretic hypotheses. The "Exponential-Time Hypothesis" (combined with results by Impagliazzo et al. [START_REF] Impagliazzo | Which problems have strongly exponential complexity[END_REF]) roughly states that 3SAT on formulas with m clauses cannot be decided in 2 o(m) time. This lower bound transfers since the constructed graph G has only 21m edges.

Theorem 1. Minimum Deletion into Disjoint Stars is NP-hard and cannot be solved in 2 o(|E(G)|) time on graphs G, even if G is restricted to maximum degree three, assuming the Exponential-Time Hypothesis.

A Heuristic for MDDS

In light of the hardness result presented in Section 3.1, we devised a heuristic that first computes a maximum-weight spanning tree T on each connected component of BG. It then transforms each T into a star packing by computing MDDS on T .

We present a dynamic programming algorithm solving MDDS on a tree T . To this end, we root T at an arbitrary vertex, and compute a dynamic programming table for vertices in a post-order traversal. Consider a set S of edges that, after removal from T , yields a collection of disjoint stars. We denote the result of this removal as T -S = (V (T ), E(T ) \ S). Each vertex x has one of three states relative to S:

D i * (x) := D i-1 * (x) + min D + (v i ), ω(xv i ) + min(D * (v i ), D -(v i ))
That is, if x is the center of a star, then the edge xv i must be in S if either v i is the center of a star or the edge between v i and its parent x is not in T -S.

D i + (x) := D i-1 + (x) + ω(xv i ) + min(D * (v i ), D -(v i ))
That is, if x is a leaf of a star centered at the parent of x, then the edge xv i must be in S.

D i -(x) := min D i-1 -(x) + ω(xv i ) + min(D * (v i ), D -(v i )), D i-1 + (x) + D * (v i )
The case of D i -(x) is a bit more subtle. Since x is not the center of a star, all but at most one edge between x and its children are in S, so if xv i is not in S then D i-1 + (x) forces all xv j to be in S, for 1 ≤ j < i. Finally, the subsolutions rooted at x are, then, given by:

D + (x) := D m + (x) D -(x) := D m -(x) D * (x) := D m * (x)
4 Quantifying Hierarchical Conficts

We applied our MDDS heuristic to homology statements on a set of prokaryotes, and on a set of eukaryotes. The solution to MDDS provides an estimate of the minimum number of positions that must be ignored so that the necessary conditions for a hierarchy, highlighted by Lemma 1, are achieved. Before applying the heuristic of Section 3.2 we cleaned the syntenic blocks according to Appendix B, and preprocessed the graphs for segmental duplications according to Appendix C.

Discordance Ratio and Distinction from Jaccard Index

Define coverage(B) of a collection of blocks B as the total number of positions covered by all segments in blocks of B. We report the hierarchical discordance ratio between collections A and B as d(A, B) = w/(coverage(A) + coverage(B)), where w is the weight of the MDDS on BG(A, B). A discordance ratio of 0.1 means that we have to ignore at least 10% of the total coverage of the blocks (in both methods) in order to have a hierarchical relationship between them. Alignathon used mafComparator to compute the straightforward Jaccard similarity index between collections of blocks. In this case, the elements of the sets in question are the pairwise alignments of positions implied by the blocks. So if a pair of positions are aligned in one class of blocks but not the other, this will contribute one to the denominator.

Consider collections A and B such that A only contains blocks with segments

from {g i [200x+1..200x+100] | 0 ≤ x < ⌊ ℓ(gi)
200 ⌋}, and B only contains blocks with segments from

{g i [200x + 101..200x + 200] | 0 ≤ x < ⌊ ℓ(gi)
200 ⌋}, for all genomes g i with length ℓ(g i ). In other words, the collections can only have blocks of length 100 that do not overlap with each other. In this case the Jaccard similarity measure will be zero no matter the length of the genomes, indicating the most severe dissimilarity, whereas the two collections are hierarchically related, showing no conflicts, and the block graph is composed only of degree zero vertices. In that sense, our comparison method is tolerant to collections that conservatively make no assertion about a region.

We consider the two measures complementary in that they capture different qualities of the overlap properties of block collections. We see in Section 4.3 instances from the Alignathon data where the Jaccard similarity is low, yet the two collections are hierarchically related, and vice versa.

Mycobacterium tuberculosis clinical isolates

For the prokaryotes, we used a set of 94 Mycobacterium tuberculosis strains [21,5] with homology statements given by the methods listed in Table 1. These sets of blocks are those produced in [START_REF] Elghraoui | Evaluating impacts of syntenic block detection strategies on rearrangement phylogeny using Mycobacterium tuberculosis isolates[END_REF], where the methods were compared to assess their impact on inferring rearrangement phylogenies. Note that all Cactus blocks used in this subsection had segments with fewer than 50 positions filtered out.

The collections of blocks for four of the methods, along with their maf2synteny counterparts are compared in Fig. 4. As expected, each collection of blocks had a very low discordance ratio with its counterpart agglomerated by maf2synteny. Further, the agglomerated blocks always have lower discordance to all the other methods, when compared to their unagglomerated counterparts. Of the unagglomerated methods, SibeliaZ is the least discordant.

There were a couple of surprises. The first is that the most discordant pairs are between the gene-based annotation method and the de novo inference methods Cactus and SibeliaZ. Contrary to the other methods, the agglomerated annotation blocks show a small improvement against Cactus, and a surprising degradation (going up from 4% to 5%) in the discordance ratio for SibeliaZ.

Table 1: Homology statement determination methods applied to the M. tuberculosis genomes

Method Description

Cactus(SNP) Cactus [START_REF] Armstrong | Progressive Cactus is a multiple-genome aligner for the thousand-genome era[END_REF] alignment guided by a ML tree based on concatenated substitutions with respect to reference strain H37Rv (NCBI accession NC 000962.3) Cactus(SibeliaZ) Cactus alignment guided by a MLWD [START_REF] Hu | MLGO: phylogeny reconstruction and ancestral inference from gene-order data[END_REF] adjacency tree computed from SibeliaZ+M2S synteny blocks.

Cactus(Mash)

Cactus alignment guided by a B(I)ONJ tree based on the genomes' Mash [START_REF] Ondov | Mash: fast genome and metagenome distance estimation using MinHash[END_REF] distance matrix.

SibeliaZ

Locally collinear blocks produced in the first step of the SibeliaZ pipeline [START_REF] Minkin | Scalable multiple whole-genome alignment and locally collinear block construction with SibeliaZ[END_REF].

Annotation

Simultaneous annotation and orthology assignment by 95% amino acid sequence identity and 95% alignment coverage.

Modifiers Description +out

Synteny blocks computed while including the outgroup strain M. canettii (NCBI accession NC 019951.1)

+M2S

Agglomerated with maf2synteny [START_REF] Kolmogorov | Chromosome assembly of large and complex genomes using multiple references[END_REF] This implies that either 1) many blocks from the Cactus method bridge between coding regions, or 2) many duplicate regions are assigned in discordant ways.

The second surprise is that the unagglomerated Cactus methods, with different guide trees, are more discordant from each other than they are with SibeliaZ. It has been reported that Cactus's sensitivity to guide trees also has implications on the downstream phylogenetic analyses [START_REF] Elghraoui | Evaluating impacts of syntenic block detection strategies on rearrangement phylogeny using Mycobacterium tuberculosis isolates[END_REF]. In Fig. 5, the checkered pattern shows that the inclusion of an outgroup affects Cactus blocks more than the choice of a guide tree. The inclusion of the outgroup strain also decreases the discordance between the Cactus blocks on different guide trees. For example, Cactus(Mash) has discordance ratios of 0.044 and 0.061 against Cactus(SNP) and Cactus(SibeliaZ), but for Cactus(Mash)+out these values are 0.022 and 0.026. Table 2 shows the discordance between a method and its version with the outgroup. Cactus is most highly affected by the inclusion of the outgroup. While SibeliaZ is somewhat affected, Annotation is barely affected. Agglomerating the blocks with maf2synteny diminishes the discordance in all cases but Annotation.

Alignathon

The Alignathon competition was created to compare "whole genome alignment" methods [START_REF] Earl | Alignathon: a competitive assessment of whole-genome alignment methods[END_REF]. Authors of WGA software were invited to submit the collections of blocks computed by their program, which were compared using the measures described in the introduction. The project fabricated two synthetic datasets that were used to evaluate the block collections, one that mimicked the properties of set of Primates, and another that mimicked the properties of a set of Mammals. We applied our MDDS heuristic to each pair of block collections. Note that we were limited to the collections available on the Alignathon downloads page, so were unable to compare to some methods, such as Mercator/Pecan [START_REF] Paten | Enredo and pecan: genome-wide mammalian consistency-based multiple alignment with paralogs[END_REF]. The results for the Primate dataset are depicted in Fig. 6 while results for the mammal dataset are depicted in Fig. 7. Being evolutionarily closely related, the Table 2: Comparison of discordance ratios between blocks computed by the same method, with and without an outgroup in the input genome set. maf2synteny usually reduces the high discordance (shown as a percentage) between the blocks. For example, Cactus(SNP) applied to the TB sets with and without the outgroup shows a very high divergence ratio, yet a much lower one after maf2synteny has been applied to the two block collections. Primates dataset mostly shows discordance ratios below 2%. While this trend is consistent with the Alignathon findings, including GenomeMatch2 (SoftBerry, Mount Kisco, NY) being relatively more discordant, there were differences with the Jaccard index reported by Alignathon. VISTA-LAGAN [START_REF] Dubchak | Multiple whole-genome alignments without a reference organism[END_REF], for instance, stands out as generally more discordant than the others, being rather dissimilar to PSAR-Align [START_REF] Kim | PSAR: measuring multiple sequence alignment reliability by probabilistic sampling[END_REF], AutoMz, and Multiz [START_REF] Blanchette | Aligning multiple genomic sequences with the threaded blockset aligner[END_REF]. EBI-MP stands out as having both the best, and the worse discordance ratios of the dataset; despite having a ratio of over 12% against progressiveMauve [START_REF] Darling | progressivemauve: multiple genome alignment with gene gain, loss and rearrangement[END_REF], it also is the only method in the set to be hierarchically related to another one (Robusta [START_REF] Notredame | Robusta: a meta-multiple genome alignment tool[END_REF]). Cactus has very low discordance with all methods except GenomeMatch2.

The simulated mammal dataset contained genomes that were separated by a larger evolutionary distance, and this was reflected in surprisingly large discordance ratios. We observe several discrepancies with the Jaccard distances reported by Alignathon ([9] -Fig. 8B). GenomeMatch3 was extremely dissimilar to all methods but Mugsy [2], yet we observe high hierarchical discordance ratios only against Robusta and Multiz. Cactus has low hierarchical discordance against all other methods, whereas it had high Jaccard distances against Mugsy, GenomeMatch3, and EBI-EPO. On the other hand, Robusta seemed to have poor comparisons for both the Jaccard and hierarchical measures.

Discussion and Conclusions

In this article we addressed the question of how to relate two collections of homology statement blocks to each other. We established a relationship between collections where we allowed overlapping parts of those collections to be hierarchically related. In the absence of these conditions, we developed a method that gives a lower bound on the number of positions that must be ignored in order for the two to be hierarchically related. The notion of being hierarchically related depends on semantics that we imposed on the blocks, which speak to the pairwise homology relationships between the constituent genomic positions appearing in the blocks. As Ghiurcuta and Moret [START_REF] Ghiurcuta | Evaluating synteny for improved comparative studies[END_REF] used "homology statements" to define their "syntenic blocks", we used "positive homology witness" pairs to limit which segments can be contained within a homology statement block; while they required every homology statement within a segment to occur in all other segments of the block, we allowed positions that do not appear in a positive witness pair in the block, as long as they do not occur in another block. We went further by associating semantic meaning to the fact that two positions appear in different blocks. This allowed us to define what it means for some (sub)collection of blocks to generalize another (sub)collection.

On the algorithmic side, we showed the Minimum Deletion into Disjoint Stars problem to be NP-Complete. Our heuristic for MDDS is based on a dynamic program that solves MDDS exactly on a tree. Future improvements will include the exploration of other algorithms with provable guarantees to the quality of their solutions. The solution to the MDDS problem gives the number of nucleotides that must be ignored so as to make the components of the block graph stars. This is a necessary condition for the two collections to be hierarchically related, but not sufficient, and thus is a lower bound on the number of nucleotides that must be ignored so as to make the two collections hierarchically related. Future work will explore ways to tighten this bound.

We studied block collections on a set of 94 Mycobacterium tuberculosis strains, built by annotation and non-annotation based means. We showed on this data that the agglomeration of blocks using maf2synteny almost always yielded collections that were less discordant. We showed surprising discordance between the gene-based annotation method and the de novo block inference methods Cactus and SibeliaZ. Cactus showed great heterogeneity, dependent on the guide tree that was used to construct the blocks.

When performing a phylogenetic analysis on the blocks, one is tempted to incorporate an outgroup for the sake of rooting the tree. We showed the inclusion of that outgroup had drastic effects on blocks, producing blocks that were less sensitive to the Cactus guide tree. This was concordant with our results from a phylogenetic study [START_REF] Elghraoui | Evaluating impacts of syntenic block detection strategies on rearrangement phylogeny using Mycobacterium tuberculosis isolates[END_REF].

We studied block collections from the Alignathon project. The simulated Primates dataset showed that EBI-MP had both the best discordance ratio, and the worst, among all pairwise comparisons, being hierarchically related to the Robusta blocks while having ratio over 0.13 with progressiveMauve. For the less closely related simulated mammalian genomes, we showed several discrepancies between the Jaccard index reported by Alignathon and our discordance ratio, the most notable one being that while Cactus had a poor Jaccard index against a few methods, it had very low hierarchical discordance with all other methods (except GenomeMatch2).

While WGA tools and syntenic block agglomeration methods have continued to be developed, the methods to compare and analyze them has lagged behind, and the definitions of syntenic blocks are usually procedural or based on colinearity. In this article we outlined constraints on homology blocks based on the homology relationships between pairs of positions in the genome. These constraints put as much importance on the ends of the blocks as it does their contents; if two genomic segments are put into different blocks, we interpret this as a statement that should only be contradicted in a generalization of the blocks. Our new measure should inform future block inference tool development, and serve as a sanity check for the practitioner studying large scale structure of sets of genomes.

Availability of Code

All of the code associated with this paper is publicly availble at the following URL: https://bitbucket.org/thekswenson/homology-evaluation.

A NP-Hardness of MDDS

Note that the notion of star and induced star coincide on bipartite graphs since, for any bipartite G, the vertices of any star-subgraph of G also form an induced star in G. Further, no collection of node-disjoint stars can contain the triangle C 3 or the path on 4 vertices P 4 as a subgraph and it can be seen that this condition is also sufficient.

Observation 1 A bipartite graph G is a collection of stars if and only if G does not contain a P 4 subgraph.

For the correctness proof, we will make two assumptions on the structure of the input formula φ, without loss of generality. First, we assume that no variable occurs in all clauses. If a variable x does occur in all clauses, then we simply add a new variable y and the singleton clause on y. Second, we assume that each clause in φ has exactly three literals. If a clause C has at most two literals, we can simply double the occurrence of any literal in C. Lemma 2. Let φ be an instance of 3SAT and let (G = (V, E), ω) be the result of applying Construction 1 to φ. Then, φ is satisfiable if and only if (G, ω) has a star packing of weight at least 12m 2 + m.

Proof. For each variable x i of φ, let us define the edge sets

T i := 0≤j<2ni {v 3j i v 3j+1 i , v 3j+2 i v 3j⊕3 i } and
F i := 0≤j<2ni {v 3j+1 i v 3j+2 i , v 3j+2 i v 3j⊕3 i }
where 3j ⊕ 3 := (3j + 3) mod 6n i . Note that any v j i has degree two in subgraph (V, T i ) if and only if j ≡ 0 mod 3 and any v j i has degree two in (V, F i ) if and only if j ≡ 2 mod 3. Further, ω(T i ) = ω(F i ) = 4mn i . We prove the two directions of the lemma separately.

⇒: Let φ be satisfiable, that is, there is a set L of literals over variables in φ such that each clause C k intersects L in at least one literal ℓ k and L contains exactly one of x i and ¬x i for all i. If ℓ k is the literal x i in clause C k = γ j i , then let e k := u k v 6j i and, if ℓ k is the literal ¬x i , then let e k := u k v 6j+2 i . Note that all e k are distinct, e k ∈ E(G) for all k. Let S clause contain e k for all clauses C k of φ and note that ω(S clause ) = m. Further, for all variables x i of φ, let S var i := T i if x i ∈ L and S var i := F i , otherwise (that is, ¬x i ∈ L). Finally, let the selected edges be S := S clause ∪ i S var i (see gray edges in Fig. 3), noting that ω(S) = 12m 2 +m. It remains to show that (V, S) does not contain a P 4 as a subgraph. Towards a contradiction, assume that (V, S) contains a P 4 p := (a, b, c, d). By construction, neither i S var i nor S clause contains a P 4 , and p must contain edges from both of these sets. Thus, p contains u k for some clause C k . Since any u k has degree one in (V, S), we can assume without loss of generality that ab = e k . By definition of e k , there are i and j such that either b = v 6j i and x i ∈ L ∩ C k or b = v 6j+2 i and ¬x i ∈ L ∩ C k . Since 6j ≡ 0 mod 3 and 6j + 2 ≡ 2 mod 3 we know that in both cases b has degree two in (V, S var i ), and both of its neighbors have degree one in (V, S var i ) and, thus, in (V, S). This contradicts (a, b, c, d) being a path in (V, S). ⇐: Let S be a maximum-weight subset of E such that ω(S) ≥ 12m 2 + m and (V, S) does not contain a P 4 as a subgraph. First, let S clause denote the set of edges of S incident with a clause node u k . Second, for each x i , let S var i denote the set of edges of S on the variable cycle corresponding to x i and note that, for each P 3 (a, b, c) in (V, S var i ), both a and c have degree one in (V, S var i ) since, otherwise, (V, S) contains a P 4 . For each i, the connected components of (V, S var i ) are paths of lengths 1, 2, or 3 and we denote the number of P 1 s, P 2 s, and P 3 s in (V, S var i ) by r i , s i and t i , respectively. By construction, each P 2 is adjacent to at most one clause vertex u k in (V, S), and since (V, S) does not contain P 4 , each P 3 is also adjacent to at most one clause vertex u k in (V, S).

Claim.

i t i = 6m.

Proof. By decomposing the 18m vertices of the variable cycles into P 3 subgraphs separated by single edges, the upper bound of 6m is attained. It suffices to show i t i ≥ 6m so, towards a contradiction, assume that i t i < 6m. Then, there is a variable x i such that t i < 2n i implying |S var i | ≤ 4n i -1 by construction. Let S ′ result from S by removing all edges incident with vertices of the variable cycle corresponding to x i and adding the edges in T i . Since x i does not occur in all clauses, we removed edges of total weight strictly less than m+(4n i -1)m = 4mn i and we added edges of total weight m|T i | = 4mn i . Since neither (V, S) nor (V, T i ) contains a P 4 , neither does (V, S ′ ), thus contradicting optimality of S. ■ Corollary 1. Each subgraph (V, S var i ) decomposes into disjoint copies of P 3 .

Corollary 2. Let v j i and v j ′ i be nodes of degree two in some subgraph (V, S var i ). Then, |j -j ′ | ≡ 0 mod 3. Corollary 3. Let u k v j i be an edge in S clause . Then, v j i has degree two in (V, S var i ).

Note that each P 3 in (V, i S var i ) has weight exactly 2m, so S contains exactly m edges of S clause . Further, by Corollary 3, all clause vertices u k have degree at most one since they are not adjacent to degree-one vertices. Together, this means that all clause vertices u k are incident to exactly one edge in S.

We now construct an assignment β and show that it satisfies φ. To this end, let β(x i ) = TRUE if and only if S contains the edge u k v 6j i for some j, k ∈ N. Note that, if S contains the edge u k v 6j i for any j, k ∈ N then, by Corollary 1, v 6j i has degree two in (V, S var i ). Then, by Corollary 2, S cannot contain the edge u k ′ v 6j+2 i for any j ′ , k ′ ∈ N. Thus, β is well-defined. It remains to show that β satisfies φ. To this end, let C k be any clause in φ, let u k z be the unique edge incident with u k in S and let x i be the variable whose variable cycle contains z. If x i occurs non-negated in C k , then z = v 6j i for some j ∈ N by construction. But then, β(x i ) = TRUE and x i satisfies C k . If x i occurs negated in C k , then z = v 6j+2 i for some j ∈ N by construction. But then, β(x i ) = FALSE and x i satisfies C k . In both cases, C k is satisfied.

⊓ ⊔

B Collections of Block that are not Clean

Many of the software we studied produced blocks that were not clean, containing blocks with overlapping segments. We removed overlapping segments by visiting pairs of blocks in an arbitrary order, removing the overlap between their overlapping segments. Although the order in which overlaps are removed can effect the final set of blocks, we made the process deterministic by visiting the pairs in a fixed order.

C Segmental Duplications

If some method makes orthology predictions that may contain multiple segments from the same genome (e.g. clusters of orthologous groups that contain paralogs from a single genome), the block graph may provide insight into how to refine the orthology groups using blocks from another method. This section outlines such a case. When a block A ∈ A contains multiple segments from multiple genomes, blocks from another set B1, B2 ∈ B could overlap in ways that create non-star graph topologies. Fig. 8 shows one such example. The blocks B1, B2 ∈ B each overlap with one of the two copies. This configuration creates the non-star topology depicted in the middle. The block A1 can easily be split into two so that the graph becomes only stars. This results in a refinement of the blocks of A, based on the blocks of B.

Blocks B1 and B2 each overlap one of the two duplicate copies of A1 in genomes g 1 and g 2 . The block A1 can be split into two blocks A1 ′ and A1 ′′ such that the collections {A1 ′ , A1 ′′ , A2, A3} and B are hierarchically related. The two connected components of BG({A1 ′ , A1 ′′ , A2, A3}, B) are both stars with vertices {A2, B2, A1 ′ } and {A1 ′′ , B1, A3}. This transformation can be generalized to vertices of higher degree, as long as the overlapping segments can be split in this way.

Fig. 1 :

 1 Fig.1: To the right, the collections of blocks A = {A1, A2, A3} and B = {B1, B2, B3} appearing in genomes g 1 and g 2 , along with the graph BG(A, B). Segments are depicted with brackets and lined up according to their positions on the chromosomes. They are labeled by their tuple (when space permits) and the block to which they belong. The configuration of positive and negative witness pairs shows that B generalizes A. Some of the genome positions are highlighted with boxes, and two such positions are connected by a solid line if they appear as a positive witness in B, and that line is bold if they are also in A. The dashed line represents one of the (many) negative homology witnesses between A1 and A2 that are negative witness pairs for A but not for B. To the left, the graph BG(A, B) appears with edges labeled by overlap length in gray. All of the connected components are stars.

2 .Fig. 3 :

 23 Fig. 3: Example of Construction 1. The clause C := (x 1 ∨x 3 ∨¬x 8 ) corresponding to the center vertex is equal to γ3 1 = γ 1 3 = γ 5 8, that is, it is the 4 th clause containing x 1 , the 2 nd clause containing x 3 and the 6 th clause containing x 8 . A truth-assignment setting x 1 to TRUE and x 3 and x 8 to FALSE corresponds to the star cover indicated by gray highlights. Note that taking the edge between v18 1 and C instead of the edge between v 32 8 and C corresponds to satisfying the clause C by x 1 instead of ¬x 8 .

1 .

 1 x is the center of a star in T -S (covered by D * (x) in the DP), 2. x has degree one in T -S and the edge between x and its parent is in T -S (covered by D + (x) in the DP), and 3. x has degree one in T -S and the edge between x and its parent is not in T -S (covered by D -(x) in the DP). Then, D * (x), D + (x), and D -(x) contain the weight of an optimal solution S x for the subtree rooted at x, for each of the three cases respectively. If x is a leaf of T , then set D * (x) := D -(x) := D + (x) := 0. Otherwise, let v 1 , v 2 , . . . , v m denote the children of x in T . We visit the children in this order and accumulate each partial subsolution, starting with D 0 * (x) := D 0 + (x) := D 0 -(x) := 0 and proceeding as follows for each 1 ≤ i ≤ m:

Fig. 4 :Fig. 5 :

 45 Fig. 4: The discordance ratio between each pair of block collections. Each of methods Cactus(SNP), Cactus(SibeliaZ), SibeliaZ, and Annotation along with their maf2synteny counterpart.

Fig. 6 :

 6 Fig. 6: Discordance ratios for simulated Primates.

Fig. 7 :

 7 Fig. 7: Discordance ratios for simulated Mammals.

Fig. 8 :

 8 Fig.8: The block A1 ∈ A has two (duplicated) segments in genomes g 1 and g 2 . The blocks B1, B2 ∈ B each overlap with one of the two copies. This configuration creates the non-star topology depicted in the middle. The block A1 can easily be split into two so that the graph becomes only stars. This results in a refinement of the blocks of A, based on the blocks of B.

  depicts two collections of homology statement blocks A = {A1, A2, A3} and B = {B1, B2, B3}. The blocks of A are A1 = {g 1 [1..12], g 2 [13..24]}, A2 = {g 1 [14..28], g 2 [26..41]}, A3 = {g 1 [53..66], g 2 [101..113]}, while the blocks of B are B1 = {g 1 [1..28], g 2 [13..41]}, B2 = {g 1 [46..68], g 2 [94..115]}, and B3 = {g 1 [34..39], g 2 [42..46], g 2 [116..121]}.

Acknowledgements

The authors would like to thank the reviwers for their helpful suggestions. AE and FV are supported by the NIAID grant R01AI105185.


Fei Zhou

Jérôme Dyment

Chunhui Tao
email: taochunhuimail@163.com

Tao Wu

Magmatism at oceanic core complexes on the ultraslow Southwest Indian Ridge: insights from near-seafloor magnetics

Oceanic core complexes (OCCs) and detachment faults (DFs) play a key role in crustal accretion at slow and ultraslow spreading centers. We investigate the effect of different magma supply at three OCCs of the Southwest Indian Ridge (SWIR) using high-resolution deep-sea bathymetric and magnetic data. The average equivalent thickness of extrusive basalt deduced from the magnetic anomalies, a proxy for the magma supply, decreases from Yuhuang (49.25°E) to Longqi (49.65°E) to Junhui (51.75°E) OCCs. Conversely, serpentinite outcrops become more abundant, the domal OCC morphology flattens as the footwall rotation (measured by the magnetization vector inclination) increases, and hydrothermal evidences become sparse. Combined with results from the amagmatic Easternmost SWIR, our study shows that the magma supply controls the character and evolution of the OCCs and DFs on the SWIR.

INTRODUCTION

As an essential constituent of slow and ultraslow spreading centers, Oceanic Core Complexes (OCCs) result in asymmetric flanks, corrugated seafloor, and mantle rocks exhumation. They have not only been found at magma-starved ridge segment ends [START_REF] Cann | Corrugated slip surfaces formed at ridge-transform intersections on the Mid-Atlantic Ridge[END_REF][START_REF] Tucholke | Megamullions and mullion structure defining oceanic metamorphic core complexes on the Mid-Atlantic Ridge[END_REF][START_REF] Blackman | Origin of extensional core complexes: Evidence from the Mid-Atlantic Ridge at Atlantis Fracture Zone[END_REF] but are widespread and account for 50% of the recent seafloor formed at the Mid-Atlantic Ridge (MAR) [START_REF] Smith | Widespread active detachment faulting and core complex formation near 13° N on the Mid-Atlantic Ridge[END_REF][START_REF] Escartín | Central role of detachment faults in accretion of slow-spreading oceanic lithosphere[END_REF]. Numerical models suggest that magmatism plays a major role in the development of the OCCs [START_REF] Buck | Modes of faulting at mid-ocean ridges[END_REF][START_REF] Tucholke | Role of melt supply in oceanic detachment faulting and formation of megamullions[END_REF][START_REF] Olive | The structure of oceanic core complexes controlled by the depth distribution of magma emplacement[END_REF]. In contrast, OCCs are still poorly known at ultraslow spreading ridges. An end-member is the Southwest Indian Ridge (SWIR) east of Melville Fracture Zone (FZ), a nearly amagmatic spreading center which exposes smooth seafloor interpreted as the remnants of DF footwalls [START_REF] Cannat | Modes of seafloor generation at a melt-poor ultraslow-spreading ridge[END_REF]. Seafloor spreading is purely tectonic there. Conversely, the SWIR exhibits a more robust magma supply further west [START_REF] Sauter | The Southwest Indian Ridge between 49°15′E and 57°E: focused accretion and magma redistribution[END_REF]. In this area, between the Indomed and Gallieni FZs (Fig. 1A), the shallower water depth and lower Mantle Bouguer Anomaly indicate a thicker crust and hotter mantle, and the higher amplitude of the Central Magnetic Anomaly a thicker extrusive basalt layer. We use high-resolution bathymetry and magnetic data to investigate three OCCs, the Yuhuang, Longqi, and Junhui OCCs, in this area, and show that the morphology, magnetic signature, and evolution of these OCCs are tightly related to the magma supply.

DATA AND METHODS

Near-bottom bathymetry and magnetic data were acquired by "Qianlong II" Autonomous Underwater Vehicle (AUV) during DY-43, 49 and 52 cruises from 2014 to 2020. The AUV collected data at an altitude of ~100 m above seafloor, with a sampling interval of 1 m along lines 400 m apart. Noise introduced by attitude variation is reduced by the method of [START_REF] Tao | Correction of tri-axial magnetometer interference caused by an autonomous underwater vehicle near-bottom platform[END_REF], and the International Geomagnetic Reference Field (IGRF) model [START_REF] Thébault | International Geomagnetic Reference Field: the 12th generation: Earth, Planets and Space[END_REF] is removed. The magnetic anomaly is gridded at 100 m×100 m interval and inverted to equivalent magnetization by using regularized minimum residual method (See details in Data Repository) assuming a magnetized layer 100 m-thick below the seafloor. Because the anomaly is measured on an uneven observation surface and for the sake of resolution, we calculate the Reduced-To-the-Pole (RTP) magnetic anomaly by forward modelling from the equivalent magnetization [START_REF] Szitkar | What causes low magnetization at basalthosted hydrothermal sites? Insights from inactive site Krasnov (MAR 16°38′N)[END_REF].

RESULTS

Geological description of the OCCs

We interpret the components of the three OCCs (breakaway, detachment surface, termination) on the shipboard and AUV bathymetry following the criteria of [START_REF] Smith | Widespread active detachment faulting and core complex formation near 13° N on the Mid-Atlantic Ridge[END_REF].

Yuhuang OCC (Fig. 1B; Suppl. Fig. DR1) strikes in ENE direction, slightly oblique to the E-W axial volcanic zone (AVZ). It extends ~6 km along the spreading direction and drops ~1.5 km in elevation from breakaway to termination. The termination separates volcanic terrains to the North and the detachment surface where serpentinite samples were dredged. Further south a chaotic volcanic area, where dredges recovered basalt, hosts a hydrothermal field with large sulfide deposits [START_REF] Tao | Seafloor hydrothermal activity and polymetallic sulfide exploration on the southwest Indian ridge[END_REF][START_REF] Liao | Bulk geochemistry, sulfur isotope characteristics of the Yuhuang-1 hydrothermal field on the ultraslow-spreading Southwest Indian Ridge[END_REF].

Longqi OCC (Fig. 1C; Suppl. Fig. DR1) extends ~15 and ~4 km along and across axis, respectively. South of the termination, the detachment surface displays corrugations in the spreading direction. Mass wasting is clearly seen on the western end of the OCC. Basalt was dredged all over the area, but serpentinite was only dredged on the detachment surface. High temperature hydrothermal site Longqi-1 was found on the hanging wall [START_REF] Tao | First active hydrothermal vents on an ultraslow-spreading center[END_REF][START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF], whereas inactive site Longqi-3 was found in the chaotic area separating the detachment surface and the breakaway. Further south, an older volcanic zone and a large DF are identified on the regional bathymetry (Fig. DR1).

Junhui OCC (Fig. 1D; Suppl. Fig. DR1), the easternmost one, extends ~9 and ~6 km along and across axis, respectively. The detachment footwall is quite flat. Mass wasting is observed in the middle of the area, and smooth seafloor is tentatively identified to the West (Fig. 1D andF). Basalt has been dredged east of the area, and serpentinite on the detachment surface only. No hydrothermal field has been reported yet.

Magnetization Inclination and Footwall Rotation

The magnetization acquired by young oceanic basalt at mid-ocean ridges within the last million year is approximately parallel to the present-day geomagnetic field. However, the magnetization vector direction can be affected by subsequent tectonics, such as DFs, as observed on the MAR at 26°N [START_REF] Szitkar | Near-seafloor magnetics reveal tectonic rotation and deep structure at the TAG (Trans-Atlantic Geotraverse) hydrothermal site[END_REF][START_REF] Szitkar | Detachment tectonics at Mid-Atlantic Ridge[END_REF]. It is not easy to estimate the footwall rotation at the observed OCCs in the absence of any specific structural or paleomagnetic measurements (e.g., [START_REF] Morris | Footwall rotation in an oceanic core complex quantified using reoriented Integrated Ocean Drilling Program core samples[END_REF]. To address this problem, we explore the full range of possible remanent magnetization directions, calculate their corresponding RTP anomalies and select the range of inclinations that make the anomalies centered on its causative source, in general basaltic features such as the linear breakaways. OCC footwall rotation axes are generally orthogonal to the spreading direction, i.e. broadly east-west for the broadly north-south spreading SWIR. Following the dipole field hypothesis, the magnetization vector trends North and its rotation around a perpendicular axis allows the rotation angle to be estimated directly from the inclination data. The inclinations that properly deskew the observed anomalies are -75°~-105° (average -90°), -75°~-120° (av. -97.5°), and -85°~-135° (av. -110°) for Yuhuang, Longqi, and Junhui OCCs, respectively. Because the inclination of the dipole field is -57° in the area, we estimate the amount of rotation as 33°±15°, 40.5°±22.5°, and 53°±25° for the Yuhuang, Longqi, and Junhui, respectively.

RTP Magnetic Anomaly

RTP anomaly over the Yuhuang OCC trends NE-SW, parallel to the structural orientation (Fig. 1B andE). The strong linear positive anomaly separates the area in three parts, from North to South (1) low and slightly negative anomaly reflecting the weak induced magnetization of recovered serpentinites; (2) strong positive linear anomaly associated with basalt, suggesting that the contrast of magnetic signature with the northern part corresponds to a lithological boundary; and (3) strong negative anomalies associated with basalt, suggesting that the transition corresponds to a polarity boundary, most likely the Brunhes-Matuyama boundary (0.78 Ma; e.g., [START_REF] Cande | Revised calibration of the geomagnetic polarity timescale for the Late Cretaceous and Cenozoic[END_REF]. The OCC therefore initiated in the late stage of the Matuyama period. RTP anomaly over the Longqi OCC is generally oriented E-W (Fig. 1C andF). In the north-west corner, around Longqi-1 hydrothermal vent, a strong positive anomaly indicates recent and abundant volcanism. A weak negative magnetic anomaly characterizes serpentinite and the detachment surface. Further south, the breakaway, a bathymetric high with basalt outcrops, is marked by three E-W lineated anomalies, from north to south a weak positive anomaly, a strong negative anomaly, and a strong positive anomaly. We interpret this succession of anomalies as a geomagnetic polarity sequence recorded by the remanent magnetization of basalt. The transitions north and south of the strong negative anomaly correspond to the polarity reversals between the Brunhes and Matuyama periods (0.78 Ma), and between the younger part of the Matuyama period and the Jaramillo normal episode (0.99 Ma), respectively. The southern side of the strong positive anomaly thus records the older end of the Jaramillo episode (1.07 Ma) and is therefore underlain by basalt although no sample is available to confirm this hypothesis. The Jaramillo anomaly further continues to the west (Fig. 1F), as interpreted by [START_REF] Wu | An intermittent detachment faulting system with a large sulfide deposit revealed by multi-scale magnetic surveys[END_REF]. The elevated area extending further south toward the termination of the older DF is also probably made of volcanics formed during the Matuyama period [START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]. Conversely, the weaker transition observed north of the weak positive anomaly corresponds to the lithological transition between serpentinites and basalt. These observations indicate that the OCC initiated during the Jaramillo episode between 0.99~1.07 Ma.

RTP anomaly over the Junhui OCC (Fig. 1D andG) shows weak negative anomalies. In the West, the inferred smooth seafloor does not display any significant signal. In the southeastern corner of the area, a negative E-W anomaly near the breakaway is interpreted as a polarity boundary, suggesting that this OCC also initiated during the Late Matuyama period. Northward, the lithological boundary between basalt and serpentinite is inferred from a weak E-W anomaly.

Modeling

Extrusive basalt is the dominant source of marine magnetic anomalies, with additional contribution of the deeper crust and/or uppermost mantle at slower spreading rates (e.g., [START_REF] Dyment | Spreading-rate-dependent magnetization of the oceanic lithosphere inferred from the anomalous skewness of marine magnetic anomalies[END_REF][START_REF] Dyment | Contribution of serpentinized ultramafics to marine magnetic anomalies at slow and intermediate spreading centres: insights from the shape of the anomalies[END_REF]. Magnetics therefore offer an indirect means to evaluate the amount of basalt and other magmatic products, and in fine to approach the magma supply for comparison among our study areas. To this end, we use forward modeling to compute the thickness of extrusive basalt that best explains the observed anomalies. However, in order to alleviate the fundamental non-uniqueness of potential field methods, we have to adopt several simplifying hypotheses to determine what we call the equivalent thickness of extrusive basalt (ETB), i.e. the thickness of extrusive basalt if our simplifying hypotheses were correct. We assume that (1) the magnetized crust is 1 km-thick beneath the seafloor -indeed deeper sources have negligible magnetic effect at 100 m altitude; (2) the magnetized crust is made of two layers: an extrusive basalt layer overlying a serpentinite layer; (3) the basalt layer bears only thermoremanent magnetization of 15A/m (See Data Repository Section D); (4) the serpentinite layer bears only induced magnetization of 2A/m (e.g., [START_REF] Oufi | Magnetic properties of variably serpentinized abyssal peridotites[END_REF]. We neglect geological features that do not have a significant magnetic effect. For instance, rubbles on the detachment surface are randomly oriented and the contribution of their magnetization to the anomaly cancels out. The deformed fault zone may bear a consistent magnetization acquired during frictional heating but its thickness is negligible.

Within the above hypotheses, we adjust the depth of the basalt-serpentinite interface, reflected by the thickness of the basalt layer, to account for the RTP anomaly. On detachment surfaces, the ETB is generally negligible, less than 100 m-thick. The polarity of the basalt remanent magnetization is set to normal (respectively reversed) where the RTP anomaly is positive (resp. negative) (See Data Repository Fig. DR3).

The modeling results indicate that the synthetic and observed anomalies are in good agreement (Fig. 2). We estimate the average amount of basalt from the ETB (See Data Repository Section E) to 77.4±10.7 m (Yuhuang), 68.5±14.2 m (Longqi), and 30.1±22.1 m (Junhui). The amount of basalt decreases gradually from Yuhuang to Longqi to Junhui OCCs. As the spreading rate is the same for the three OCCs, these values are representative of the magma supply. The Yuhuang OCC is more magmatic than the Longqi OCC, and much more magmatic than the Junhui OCC.

The ETB varies in average (on sections) from 16% (Yuhuang and Longqi) to 10% (Junhui) of the 500 m-thick basalt layer of standard oceanic crust, suggesting that this proportion of the crustal accretion is magmatic in the study area. This estimate only reflects the crustal accretion on the SWIR southern flank, as no clear OCC appears on the shipboard multibeam bathymetry of the northern flank for the considered period. In absence of high resolution bathymetric and magnetic data, no inference can be drawn for this flank, and therefore the whole area, for the last ~1.5 Ma.

DISCUSSION

Our observations and result suggest that the geology and magnetic signature of the OCCs are causally linked to the magma supply. In decreasing order of the magma supply, from Yuhuang to Longqi to Junhui OCCs, we observe (1) gradually flattening footwall surfaces as reflected by the increasing footwall rotation angles, (2) more abundant serpentinite outcrops, and (3) less evidences of hydrothermalism, whereas (4) the RTP anomalies related to magnetic polarity become weaker, less linear, and cover less area. For even lower magma supply, previous studies of the easternmost SWIR in a nearly amagmatic environment [START_REF] Cannat | Modes of seafloor generation at a melt-poor ultraslow-spreading ridge[END_REF][START_REF] Cannat | On spreading modes and magma supply at slow and ultraslow mid-ocean ridges[END_REF][START_REF] Sauter | Continuous exhumation of mantle-derived rocks at the Southwest Indian Ridge for 11 million years[END_REF] showed smooth footwall surfaces made of serpentinite on both flanks, with consistently weak magnetic anomalies. Figure 3 presents four schematic sections in decreasing order of magmatism. The more magmatic Yuhuang OCC (Fig. 3A) is probably underlain by strong gabbroic cores, resulting from the capture of gabbroic bodies by the footwall near the Lithosphere-Asthenosphere boundary, as reflected by its broad dome morphology and smaller footwall rotation [START_REF] Macleod | Life cycle of oceanic core complexes[END_REF]. At the less magmatic Longqi OCC (Fig. 3B), the footwall may include less gabbro bodies resulting in a less-marked dome structure [START_REF] Ciazela | Mantle rock exposures at oceanic core complexes along mid-ocean ridges[END_REF], although hot gabbro bodies have not been detected so far by seismic methods in this area [START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]. Furthermore, the observation of two OCCs and an intermediate volcanic zone supports the alternance of detachment and magmatic spreading. The even less magmatic Junhui OCC (Fig. 3C) shows a flatter and longer serpentinite outcrop, suggesting that the footwall only captures rare gabbroic bodies. For amagmatic spreading as the easternmost SWIR (Fig. 3D), the footwall of the DF is not even defined as an OCC sensu stricto because it does not display the typical dome morphology due to its exclusive serpentinite composition. There, DFs are purely tectonic, with the previously active DF being terminated by the initiation, across its footwall, of a new DF with opposite slope [START_REF] Sauter | Continuous exhumation of mantle-derived rocks at the Southwest Indian Ridge for 11 million years[END_REF].

Our results suggest that magmatism and tectonics are controlling in various degrees the Yuhuang, Longqi, and Junhui OCCs, which all initiated in response to a decreasing magma supply at the ridge axis and will probably end as magmatism gradually resumes. The magma supply affects characters of the OCCs such as morphology, serpentinite outcrops and underlaying gabbro bodies, hosted hydrothermal vents, and the magnetic signature. Magmatism in this area results in the initiation and termination of detachments separated by magmatic spreading periods, in alternance. Conversely, amagmatic spreading results in the purely tectonic exhumation of serpentinite on alternating DFs, as observed in the Easternmost SWIR. The variability of OCCs and DFs along the ultraslow spreading centers may reflect the heterogeneity of the along-axis magma supply and therefore offers a proxy to this parameter.

Figure 2. Forward modelling of the RTP magnetic anomalies on the three investigated areas. Each panel A-C and G-I shows, on top, a section across the structure resulting from forward modelling where green is serpentinite, red (blue) basalt with normal (reversed) polarity; and, on bottom, modelled (red) and observed (blue) RTP anomalies. The location of these section is indicated as black lines labelled 1 (panels A-C) and 2 (panels G-I) on maps in the middle column. These maps display the modelled RTP anomaly (color) superimposed on the bathymetry (shadows). See text for details on the forward modelling. Figure 3. Synthetic cross-sections presenting OCCs and DFs from more to less magmatism as suggested by the averaged equivalent thickness of extrusive basalt. A), B), C) display interpretations of the Yuhuang, Longqi and Junhui OCCs. D) displays the interpretation of the amagmatic Easternmost SWIR after [START_REF] Cannat | On spreading modes and magma supply at slow and ultraslow mid-ocean ridges[END_REF].

ACKNOWLEDGMENTS

We thank all the members of DY-43,49, 52 cruises and the AUV team. The data collection was supported by 2018YFC0309901 and 2018YFC0309903 and by the COMRA Major Project DY135-S1-01-01. Fei Zhou is funded by China Scholarship Council 201906170047. We thank Editor J. Dickens, A. Morris and four anonymous reviewers for their contribution to improve the paper. Correspondence on methodology and results of the magnetic study should be sent to JD, and on data availability and the regional project to CT. This is IPGP contribution 4256.