Scaling global tractography with ExaTract: from large cohorts to ultra-high-resolution datasets

Scaling global tractography with ExaTract: from large cohorts to ultra-high-resolution datasets Simon Legeaya,b , Felix Matuschkec,d , Bastien Herlina,e , Jan André Reuterc,d , Christophe Destrieuxf,g , Markus Axerc,d,h , Ivy Uszynskia,b , Cyril Poupona,b

arXiv:2611.03419v1 [eess.IV] 2 Nov 2026

a Université Paris-Saclay, CEA, CNRS, Neurospin, UMR 9027 BAOBAB, GAIA Laboratory, Centre d’études

Paris-Saclay, bat 135, Gif-sur-Yvette, 91191, France b AIDAS, Gif-sur-Yvette, 91191, France c Institute of Neuroscience and Medicine (INM-1), Forschungszentrum Jülich GmbH, Wilhelm-Johnen-Straße, Jülich, 52428, Germany d AIDAS, Jülich, 52428, Germany e Sorbonne University, Frontlab, CNRS UMR 7225, INSERM U1127, Hôpital Pitié-Salpêtrière, 47 bd de l’Hôpital, Paris, 75013, France f University of Tours, INSERM, Imaging Brain & Neuropsychiatry iBraiN U1253, 10 boulevard Tonnellé, Tours, 37032, France g CHRU de Tours, 2 boulevard Tonnellé, Tours, 37044, France h School of Mathematics and Natural Sciences, University of Wuppertal, Gaussstraße, Wuppertal, 42119, Germany

Abstract Global tractography offers a theoretical framework for reconstructing white matter connectivity from diffusion MRI by jointly inferring all streamlines and their consistency with the diffusion signal and underlying anatomical priors. However, its sequential nature and log-linear complexity make it computationally prohibitive at the scales demanded by neuroimaging, from large population cohorts such as the Human Connectome Project (HCP) to ex vivo ultra-high-field acquisitions at mesoscopic resolution. In this work, we present ExaTract, a parallelized spin-glass tractography framework designed for high-performance computing (HPC) architectures. Implemented as an open-source MPI+Kokkos application, ExaTract partitions the tractography domain into spatially overlapping regions processed across distributed memory nodes, with regular synchronization of boundary modifications to preserve global convergence guarantees. We validate ExaTract on both in vivo and ex vivo diffusion MRI data: the complete HCP dataset demonstrates population-scale throughput, while the high-resolution Chenonceau dataset provides the first mesoscopic structural connectome derived from global tractography.

1. Introduction Tractography is a powerful technique for reconstructing the white matter pathways in the brain using diffusion MRI data. It has been widely used in neuroscience research and clinical applications to study brain connectivity and characterize various neurological disorders. Contrary to streamline tractography, which proceeds locally with the inference of axonal fibers independently, global tractography methods aim to solve the ill-posed problem of reconstructing the white matter connectivity from the local diffusion signals by considering the simultaneous inference of all streamlines (Poupon et al., 2000; Cointepas et al., 2002; Kreher et al., 2008; Reisert et al., 2009). As emphasized by Mangin et al. (2013), this interest in global approaches

has since become increasingly reflected with the large adoption of semi-global streamline filtering methods such as SIFT2 (Smith et al., 2015) and COMMIT (Daducci et al., 2015). Global tractography, rather than performing independent streamline reconstruction and filtering as separate steps, combines both within a single inference framework, yielding tractograms that are consistent with both the measured diffusion process and anatomical priors on a global scale. Spin-glass tractography is a concrete implementation of global tractography (Poupon et al., 2000). Commonly based on a simulated annealing and Metropolis-Hastings algorithm, configurations of spin-glasses, corresponding to small portions of streamlines, are drawn from a proposal distribution and accepted or rejected based on a global energy function. The approach allows for the incorporation of global information and constraints, as the bending at white/gray matter (WM/GM) interfaces (Teillac et al., 2017) or the optimization of the diffusion model during spherical deconvolution (Kreher et al., 2008; Reisert et al., 2014; Konopleva et al., 2018). In addition, the flexibility of the energy function in global tractography allows for the adaptation of the algorithm to any data modalities, making it a powerful tool for reconstructing brain connectivity at multiple scales (Axer and Amunts, 2022). Despite these conceptual advantages, global tractography is computationally challenging. Its sequential nature creates memory and time bottlenecks that become prohibitive at the resolutions and dataset sizes encountered in practice. Large in vivo cohorts such as the Human Connectome Project (HCP) (Van Essen et al., 2013) and the UK Biobank (Sudlow et al., 2015) provide millimeter resolution data across thousands of subjects, requiring efficient processing at the population level. At the other extreme, ex vivo ultra-high field diffusion MRI at the mesoscopic scale, and microscopy technique modalities, including 3D polarized light imaging (3D-PLI) (Axer et al., 2011), polarization-sensitive optical coherence tomography (Wang et al., 2018), and scattered light imaging (Menzel et al., 2021), enable the exploration of anatomical details of the WM pathways while disentangling the common biases of in vivo diffusion MRI, such as partial volume effects, gyral bias, and fiber crossings. The resulting Big Data, combined with the log-linear complexity of the global tractography algorithm, requires its implementation to scale up in order to tackle the aforementioned bottlenecks. In this work, we introduce ExaTract, a parallelized implementation of the global spin-glassbased tractography technique designed to leverage high-performance computing (HPC) resources. ExaTract is implemented as an MPI+Kokkos application optimized for distributed memory architectures, enabling tractography at scales ranging from large population cohorts to a deep phenotyping approach at ultra-high resolution. The parallelization strategy relies on a spatial decomposition: the tractography domain is partitioned into sub-regions processed semi-independently across computing nodes, with overlapping ghost-cell regions and a communication protocol that synchronizes modifications across partition boundaries. This design is thought to preserve the global asymptotic consistency essential to the system convergence while enabling efficient parallel execution. We demonstrate the accuracy and the scalability of ExaTract on datasets representative of both in vivo population imaging and post-mortem exploration, covering the wide and deep phenotyping Big Data challenges. 1.1. Related works A recent work from Kruper and Rokem (2026) has proposed a survey of massively parallelized brain tractography frameworks pointing out the key information to be taken into account to design a suitable framework for HPC tractography. The main challenge in parallelizing the global tractography algorithm is the inherent sequential nature of its underlying Monte Carlo Markov chain. Further adaptations and assumptions were required to implement a parallelized 2

architecture while maintaining acceptable efficiency. To date, a few models have been proposed to parallelize global tractography on shared-memory architectures. The parallelization opportunity comes from the observation that, among the possible modifications, the energy variations and the corresponding reconfigurations of spin-glasses remain local, even though the energy terms remain global. Philips et al. (2013) and Wu et al. (2016) exploited this assumption by adopting an explicit parallelized global tractography scheme with two different approaches. Philips et al. (2013) relied on a parcellation scheme to define disjoint, distant areas of authorized modification, each area handled by an independent thread. The distance to other areas prevents any corrupting interaction. At regular user-defined intervals, all threads stop their work and switch to the next area, introducing a cyclic scheme. Wu et al. (2016) used a dynamic assignment of these areas. At regular intervals, the areas are dynamically recomputed. While these two techniques prevent both alteration of the cost function estimate and race conditions when two threads try to write a variable at the same time, they show limited scalability, as the effective area of processing for the global tractography is divided by more than the number of parallel workers. Konopleva et al. (2018) proposed an alternative strategy based on locking, which allows threads to safely access and modify the shared resource asynchronously. A computational challenge is the overhead introduced by mutex strategies, which increases significantly as the number of threads increases. In addition, it requires careful design to avoid deadlocks, in which two or more threads wait indefinitely for each other to release locks. This constraint practically limits the acceleration to a few threads at most. To further leverage distributed-memory architectures, an adaptation of the global tractography algorithm is necessary. The very high dimensionality of the problem, comprising a few million spin-glasses (each with 6 variables), makes it a candidate for an asynchronous strategy based on spatial decomposition, where the tractography domain is shared between the computing units. While never tried in global tractography, the massive parallelization of simulated annealing has already been demonstrated for certain applications (Jayaraman and Darema, 1988; Greening, 1990). 2. Materials and Methods 2.1. Global tractography framework Spin-glass tractography simultaneously infers the entire tractogram through an optimization process constrained by the local orientation distributions to overcome the local perturbations. First introduced by Poupon et al. (2000) and further formalized by Cointepas et al. (2002) and Mangin et al. (2002), particle-based global tractography involves particles, called spin-glasses, that represent small portions of streamlines. Each spin-glass extremity can connect with another spin-glass extremity to form chains of spin-glasses, ultimately leading to complete streamlines. 2.1.1. System and energy definition A spin-glass Xi is defined by its position xi , its orientation ni , and its half-length l, here defined as the same for all the spin-glasses. Its extremities are designated by Xiαi where αi = +1 for the one along ni or αi = −1 for the other along −ni . A generic extremity will then be denoted Xi± . The entire system M is hence composed of

3

n o   X : X0 ...XN−1     M :  n o    C : (X ± , X ± )...(X ± , X ± ) i j k l

the spin-glasses, (1) and their connections.

In the following, D denotes the data or the local orientation map, whichever local modeling is used. Historically, spin-glass-based global tractography relies on optimizing two energy terms, internal and external, such as the total energy of the system E(M , D) is given by: E(M , D) = αEint (M ) + βEext (M , D),

(2)

where α and β represent an empirical weighting of the energies. The internal energy tends to regularize connections between spin-glasses with respect to the expected axonal fiber anatomy, by optimizing the distance and incident angles between two connected spin-glasses. The energy term between two extremities, Ucon defined in Reisert et al. (2009) is commonly used as a more permissive and exploratory energy: Ucon (X1α1 , X2α2 ) =

∥x1 + α1 ln1 − x̄∥2 + ∥x2 + α2 ln2 − x̄∥2 − L, (l/2)2

(3)

where x̄ is the midpoint between the two spin-glasses and L is the connection likelihood. A large L value strongly encourages the connection, while a small one promotes only connections between very close spin-glasses. The total internal energy of the system is the sum of all the energy potentials corresponding to single connections: X Eint (M ) = Ucon (Xi± , X ±j ). (4) (Xi± ,X ±j )∈C

The external energy constraints the orientation of the spin glasses to align with the anatomical profile stemming from the measured data (dMRI, PLI, or ps-OCT...). In ExaTract, nongenerative approaches were chosen, for which the alignment of any spin glass can be expressed as a function of the log-probability of its orientation stemming from the corresponding orientation distribution function (ODF), as defined in Poupon et al. (2000): X Eext (M , D) = − log (p(Xi )), (5) Xi ∈X

where p(Xi ) is the probability to have an orientation ni in the ODF at position xi . A further regularization term µ | X | is added to the energy in order to control the number of spin-glasses and avoid any divergence of the optimization process. 2.1.2. Optimization process Global tractography generally starts with a pseudo-random spin-glass configuration of the system. ninit spin-glasses are randomly spread across the volume and oriented according to a Gibbs distribution sampler built upon the probabilities stemming from the local ODF, further sharpened or smoothed using a temperature T gibbs . Through a reversible-jump Monte Carlo Markov chain (RJMCMC) process (Green, 1995), infinitesimal modifications are sequentially proposed to the system. Unlike traditional MCMC algorithms (Neal, 1993), RJMCMC allows 4

the system to change its dimensions, meaning that a proposal can switch from one parameter space to another of a different dimension. This property makes RJMCMC particularly suitable for global tractography (Reisert et al., 2009), as the optimal number of particles needed to fit the data best is unknown in advance. Propositions from the system configuration M to M ′ are drawn by randomly selecting a spin-glass and applying modifications, including: • motion: the spin-glass is randomly moved and reoriented, as in Konopleva et al. (2018). These random adjustments are selected from a normal distribution around the current position and orientation of the spin-glass, with a standard deviation σ = 0.05. Low temperature allows local exploration of the configuration space, enabling the system to refine the alignment of spin-glasses with the local orientation distribution while allowing some degree of freedom to escape local minima. • connection: all neighboring spin-glasses with a free extremity within a sphere of defined radius r search are considered for possible connection (see Figure 1). In addition, a noconnection possibility is also tested. This draws a distribution in which the connection proposal is chosen, l rsearch

Figure 1: Connection of a spin-glass. When proposing a connection, a random extremity of a random spin-glass, possibly already connected, is drawn. All possible connections with other spin-glass extremities found in a sphere of radius r search are tested. In addition, a no-connection possibility is also evaluated. The proposal is selected among the set of tested connections with a Gibbs sampler of temperature T gibbs based on their energy.

• creation: a new spin-glass is randomly created within the tractography domain. This spin-glass orientation is chosen probabilistically using a Gibbs sampler similar to the one used for the initialization of the spin-glass map, ensuring that the proposed configuration is likely to be accepted while also maintaining enough randomness in the creation of new spin-glasses to explore a wide range of plausible configurations. To maximize the proposition likelihood, a connection is jointly proposed for the new spin-glass, which then has to be validated depending on the acceptance ratio R. • deletion: the spin-glass is deleted, releasing its external energy and its connections. If p prop (M ′ |M ) designates the proposal distribution to switch from M to M ′ , the modification of the configuration is proposed with an acceptance ratio of R, according to Bayes’ rule and under the assumption of a Gibbs distribution of the system at temperature T: R=

p prop (M |M ′ ) − ∆E e T . p prop (M ′ |M ) 5

(6)

The acceptance ratio only depends on the proposal distribution and the variation of energy ∆E between M ′ and M . In the underlying Metropolis-Hastings optimization process, which is the sampling algorithm typically used to draw a Markov chain, a modification is accepted if R ≥ 1. If R < 1, the modification is accepted with a probability of R. This process is repeated a large number of times until a stopping criterion is met, following a temperature-cooling schedule. The performed modifications tend to minimize the system energy E(M , D), leading the system configuration to converge towards an optimal tractogram that is assumed to provide a plausible configuration both from the point of view of the anatomy and the data. The final tractogram is then extracted from the optimal configuration of spin glasses by following the connections between them. 2.2. Parallelization strategy of ExaTract Each modification brought to a spin-glass Xi , initially connected to spin-glasses X j and Xk and potentially newly connected to spin-glasses X j′ and Xk′ , can be rewritten locally under the following general form: ∆E = Ucon (X ′ +i , X ±j′ ) − Ucon (Xi+ , X ±j ) + Ucon (X ′ −i , Xk±′ ) − Ucon (Xi− , Xk± )
− log p(X ′ i ) + log (p(Xi ))

(7)

+ γµ , where γ is 0 in case of motion or connection modification, 1 in case of creation, and -1 in case of deletion. This local formulation of the energy variation is a key property that enables parallelization of the global tractography algorithm. 2.2.1. Domain decomposition To fit with current parallel architectures, the tractography domain is divided into N parcel precomputed static sub-domains using the METIS library (Karypis and Kumar, 1997). The binary mask corresponding to the tractography domain can be represented as a graph, in which each voxel corresponds to a vertex and each edge represents a neighbor relationship between two adjacent voxels. The multilevel k-way partitioning algorithm (Karypis and Kumar, 1998) provided by the METIS toolbox decomposes this graph, taking into account the spatial locality of the data. The partitioning process aims to create sub-domains that are as independent as possible by following the morphology of the tractography domain while minimizing the number of edge cuts between sub-domains. Each partition being assigned to a specific shared-memory computing unit, the computed domain decomposition directly contributes to minimizing communication costs between nodes during spin-glass optimization. A strong constraint in the partitioning process is to ensure that the sub-domains are sufficiently large to allow for efficient parallel processing while also being small enough to fit within the memory of each node, which further adds to the complexity of the parallelization of global tractography. To ensure system continuity across partition boundaries, an overlap must be defined between adjacent sub-domains. This overlap allows information exchange between nodes and ensures that modifications made in one sub-domain can influence neighboring sub-domains, which is crucial to preserve the continuity of the underlying Markov field. A set of ghost cells (Kjolstad et al., 2010) has to be defined within each overlapping area and has to be updated at regular intervals to 6

reflect the changes made in the neighboring sub-domains. The width of the overlap is determined by the spin-glass length and the maximum interaction radius of spin-glasses, ensuring that all relevant interactions are captured while minimizing communication overhead (see Figure 2). A

B

Partition i

C

Partition i

l rsearch

Spin-glass interaction domain

Overlap of 1 voxel size

Overlap of 2 voxels size

Figure 2: 2-D examples of overlap depending on l, r search (A) and the voxel size. When l + r search is lower than the voxel size (B), the overlap pattern includes the 26 directly adjacent voxels. When l + r search is greater than the voxel size (C), the external overlap of the partition i needs to be at least two voxels wide. The same logic applies symmetrically on the other side of the boundary, so the minimum size of the full overlap is 2 voxels.

Gathering modifications in the overlap areas at regular intervals is a necessary assumption for adapting the underlying RJMCMC process to a shared memory environment, based on the observation that a single modification affects only a local area of the system, even if the energy terms are global. By regularly updating ghost cells at a certain frequency, we can ensure that modifications made in one sub-domain are reflected in neighboring sub-domains, ensuring coherent convergence across overlap areas and, consequently, across the entire domain. The frequency of communication across the node is then defined using an additional methodological parameter, ω, representing the workload per iteration, and corresponding to the number of proposals performed before synchronizing the ghost cells. The workload ω is user-defined in order to adapt to the size of the subdomain and ensure a balanced advancement of the system across the entire domain. ω is directly related to the communication frequency, as a higher workload corresponds to a lower communication frequency. This approach allows us to maintain the system’s global coordination while enabling efficient parallel processing of the tractography algorithm. 2.2.2. Multithreading In addition to domain decomposition, each sub-domain assigned to a node can leverage the node’s multithreading capabilities to perform concurrent execution of the tractography algorithm within that sub-domain. While previous implementations based their parallelization on scheduling and synchronization of the modifications (Philips et al., 2013; Wu et al., 2016), or mutual exclusion methods (Konopleva et al., 2018), we introduce an additional, but less restrictive assumption to allow threads to make changes in parallel without incurring any additional overhead. Threads independently draw modifications and apply them concurrently to the spin-glass maps. This assumption is expected to be less stringent than the previous one done in domain decomposition. 7

In computer science, modifying the same memory location at the same time is called a race condition. It may lead to undefined behavior, as the memory location may not hold the intended value. In ExaTract, we avoid the potential fatal consequences of race conditions by a posteriori corrections: before communications, the spin-glass map configurations are checked to correct for abnormalities. While an additional linear-complexity step is required to implement it, the added value of thread asynchronism is relevant to minimize computing time. A local and reduced copy of the system M is developed in ExaTract to mitigate race condition occurrence. The different modification proposals are evaluated on parallel branches, each operating on a private and temporary copy of the current spin-glass configuration around the spin-glass being modified. As illustrated in Figure 3, this approach minimizes the duration of concurrent access tc to shared data while fully preserving parallel execution. A model for the distribution of conflict occurrences is presented in Appendix A. Spin-glass map (or master) local copy of the involved spin-glasses local copy of the involved spin-glasses

Computation

on core 2 Computation on core 1

Have initial spin glasses changed on master ?

?

Have initial spin glasses changed on master ?

tc

no

merging yes

abandon

no

yes

? abandon

tc

merging

Figure 3: Local copy system for parallel computing of modifications. This system includes (a) a minimal copy of the involved spin-glass, (b) parallel computation on the different threads, and (c) a quick check of whether the original spinglass has changed. While allowing for parallel processing, a concurrent situation can only occur during the time lapse tc .

2.2.3. Conflict resolution The same correction mechanism is applied in ghost cells, where similar conflicts can occur when merging modifications stemming from two adjacent nodes. The conflict resolution follows a deterministic, symmetrical strategy to ensure the same correction is applied to all nodes sharing the ghost cell (for more details, see Appendix B). This approach enables rapid recovery from conflicts by efficiently restoring the system to a consistent state. Taking these anomalies into account at runtime allows us to preserve the asymptotic convergence towards an optimal configuration of the spin-glasses. However, conflicts should be minimized as much as possible, as they break the Markovian framework. This 8

minimization is possible by reducing ω, the number of operations performed per iteration, even though it also plays in favor of an increase in communication events. It raises the question of the optimal balance between communication and conflict occurrence, which is crucial for the computational efficiency of the parallelization strategy but also for the convergence of the system, as conflicts are expected to slow down the convergence. Kokkos kernels

Domain decomposition

Equivalent virtual topology on a cluster of 4 nodes

Optimisation

Chronogram

Communication Conflict correction

MPI iteration

Spin-glass container expansion

…

0

…

1

…

2

save on disk

communication

…

3

time

Orientation distribution function map

Spin-glass modifications: - connection - birth/deletion - motion -

Initialization step

Runtime transitory steps

Final step

Figure 4: Chronogram of ExaTract execution. A brain, along with its ODF map, is divided into four partitions, each assigned to a node of a supercomputer. The four MPI processes run global tractography in a synchronized way, performing a user-defined number of operations defined by ω. An MPI iteration comprises (1) an optimization kernel that effectively applies the global tractography on its partition, (2) communication steps that exchange the modifications that occurred in overlapping areas, (3) a conflict correction step, and (4) the spin-glass map expansion.

2.3. Implementation and monitoring ExaTract is fully implemented in the C++ language using object-oriented programming. It employs data structures, IO readers and writers, and dMRI processing objects, including diffusion ODFs and subsequent methods, as part of the Ginkgo C++ toolkit1 . MPI communication is managed with Isend and Irecv, which are asynchronous, nonblocking methods that allow communication to overlap with other computation. Multithreading is implemented using Kokkos (Carter Edwards et al., 2014), an open-source portability library. Kokkos abstracts parallel execution through two core concepts: kernels, which encapsulate parallel loops, and views, which manage memory allocations across architectures. ExaTract required several Kokkos kernels to execute its parallel workload: spin-glass map initialization, global tractography optimization, a posteriori correction, merging communicated 1 https://framagit.org/cpoupon/gkg

9

modifications, and a set of kernels for dynamic expansion of the spin-glass map. Because the total number of spin-glass segments per voxel is not known at the beginning, ExaTract supports dynamic allocation of new spin-glasses within the Kokkos view representing the spin-glass map. It is designed to favor cache efficiency when iterating over neighboring spin-glasses, which is especially necessary for high-resolution datasets. Algorithm 1 and Figure 4 summarize the chronological execution of these kernels. ExaTract comes with additional features, including a logging system for real-time monitoring of algorithm convergence and error occurrence analysis. Regular save checkpoints are done as ExaTract might need several consecutive jobs to reach convergence when targeting highresolution datasets. In addition, the backup ensures resilience in case of an unexpected job failure. Algorithm 1 Simplified pseudocode of the ExaTract process 1: for it = 0 → itmax do 2: for draw = 0 → ω ∗ Nvoxel /Nthreads do 3: draw a random spin-glass 4: draw a modification proposal 5: if p ≤ R then 6: apply modification 7: for spin-glass in local map do 8: if spin-glass is abnormal then 9: correct spin-glass 10: for spin-glass in local overlap do 11: if spin-glass has changed then 12: add spin-glass to message 13: for all neighbor ranks i do 14: send message 15: receive messagereceived 16: for spin-glass in messagereceived do 17: apply modification spin-glass 18: for voxel in local map do 19: expansion of voxel if necessary 20: T = T + δT

▷ per MPI process ▷ per thread

▷ parallelized over threads

▷ parallelized over threads

▷ parallelized over threads

2.4. Experimental protocols To assess the consistency and the scalability, ExaTract was tested on both synthetic and realcase datasets of varying sizes (see Table 1). ExaTract has been designed to be compatible with any kind of local model. In the frame of this application, we used the analytical Q-ball model (Descoteaux et al., 2007) computed with a Laplace-Beltrami regularization factor of 0.006. All ODFs were sampled using 10,000 uniformly distributed orientations over the unit sphere, corresponding to a solid angular resolution of 2.0 degrees. For all experiments, creation, deletion, motion, and connection probabilities ratios were set to 0.25, 0.02, 0.25, and 0.15, respectively. The connection likelihood L was set to 0.5, and the regularization of spin-glass count µ was set to 0.1. The internal and external energy weighting α/β was set to 5.0/4.0. The spin-glass length was 10

typically set slightly below the resolution in order to minimize overlap size. These parameters were either derived from the literature (Reisert et al., 2009; Konopleva et al., 2018) or empirically defined based on preliminary experiments to ensure convergence of the system while allowing for exploration of the configuration space. 2.4.1. Synthetic data The DiSCo MRI phantom (Rafael-Patino, Girard et al., 2021) was used as a synthetic dataset to evaluate the performance of ExaTract with different numbers of threads and MPI processes on the same dataset. Simulated data of a shell with a b-value of 3000 s.mm−2 was used, and the tractography domain was defined as the whole phantom. The following parameters were used for ExaTract: 10 initial spin-glasses per voxel, spin-glass length and connection radius of l = r search = 0.99 mm, 500 MPI iterations with ω = 10 operations per voxel and per MPI iteration, an exponential decreasing cooling schedule from T init = 0.4K to T f inal = 0.03K. 2.4.2. In vivo cohort The 1065 subjects of the HCP database (Van Essen et al., 2013) were processed with ExaTract using the following parameters: 10 initial spin-glasses per voxel, spin-glass length and connection radius of l = r search = 1.249 mm, 14000 MPI iterations with ω = 1 operation per voxel and per MPI iteration, exponential decreasing cooling schedule from T init = 0.4K to T f inal = 0.01K. The preprocessing comprised the application of the FreeSurfer2 pipeline to compute masks and the Ginkgo pipeline3 to compute a map of ODFs as described in Herlin et al. (2024). 2.4.3. Post-mortem high-resolution datasets The Chenonceau dataset (Legeay et al., 2026) consists of an entire cerebrum scanned with a preclinical MRI system at 11.7T (Bruker Biospin, Ettlingen), at an isotropic resolution of 200 µm according to a multiple shell sampling of the q-space including 25 directions at b = 1000s.mm−2 , 60 at b = 4500s.mm−2 and 90 at b = 8000s.mm−2 . The tractography domain was defined as the entire brain; the domain was split into 800 parcels of 144414 ± 7 voxels. The ODF map was computed using the 90 directions of the highest b-value shell. An additional dataset of a macaque brain (Macaca Fascicularis species) scanned at 250 µm isotropic resolution with the same diffusion MRI protocol (Poupon et al., 2025) was also used to evaluate the scalability of ExaTract. Dataset Type Isotropic resolution Dataset size (voxels)

DiSCo phantom 1 mm 15040

HCP subject in vivo human 1.25 mm 5 − 8x105

Macaque ex vivo specimen 250 µm 3.26x106

Chenonceau post-mortem human specimen 200 µm 1.15x108

Table 1: Sizes of the different use-case datasets.

2.4.4. Hardware For benchmarking on the DiSCo dataset, ExaTract was run using a desktop machine equipped with 8 bi-threaded CPU cores (2.40GHz, cache 16Mo). For all other datasets, the Irène-Rome 2 https://surfer.nmr.mgh.harvard.edu/ 3 https://framagit.org/cpoupon/gkg/

11

partition of the Joliot-Curie supercomputer of the TGCC (GENCI, Très Grand Centre de Calcul, CEA, Bruyères-le-Châtel, France) was used, mobilizing up to 200 Non-Uniform Memory Access (NUMA) nodes of 64 bi-process cores (256 Go RAM per node, 2.60GHz) to compute the global tractography. Due to the non-uniform memory access in this architecture, the number of MPI processes per node was typically set to 8, favoring local memory of the multi-threaded processes, as further demonstrated in section 3.1.2. 3. Results 3.1. Parallelization assumption validation The design of ExaTract relies on two main assumptions to allow for efficient parallelization: (1) the modifications performed in a voxel only impact a local area of the system so that domain decomposition is effective, and (2) the occurrence of conflicts in multithreading and in overlapping regions remains negligible in the convergence of the system. To validate these assumptions, we first estimated the mean conflict occurrence in both cases and then compared it with the observed occurrences across the different benchmarks. Based on equations A.2 and A.3, we estimated the upper limit of error occurrences from the two parallelization layers (Figure 5). The estimation of communication conflict occurrences is based on the probability that a proposed modification will affect another modified spin-glass in the overlapping area of another node. The estimation of multithreading conflicts is based on the probability that two threads draw and modify the same spin-glass, which is directly related to the number of threads and the number of operations per MPI iteration. Both estimations 2 are weighted by the squared mean acceptance ratio R , since a conflict can occur only if both proposed modifications are accepted.

12

2.5

1e 6 Dynamic =1 Dynamic =5 Dynamic =10

2.0 1.5

10 4

Error

10 5

1.0 Communication =1 Dynamic =1 Communication =5 Dynamic =5 Communication =10 Dynamic =10

10 6 0.5 10 7

0.0 0

20 40 60 Number of threads Nt

0

20 40 60 80 100 Spin-glass count per voxel n

Figure 5: Estimation of error occurrences per voxel evaluated over multiple ω values. On the left, the relationship between the multithreading conflict occurrences and the number of threads Nt shows the convergence of equation A.3 to an upper bound weighted by ω at high thread counts. On the right, dependencies of the two conflict sources on the number of spin-glasses per voxel show the preponderance of communication conflicts over multithreading conflicts (or dynamic conflicts) (Nt was fixed to 16). For the numerical applications, R was fixed to 0.1, tc to 0.001 and poverlap to 0.05. Communication conflicts occurring are counted per voxel in overlapping areas.

Numerical applications of the equations established in Appendix A, exposed in Figure 5, show that the errors stemming from domain decomposition are more likely to occur than the ones stemming from multithreading. A striking result is that multithreading conflicts converge to an upper bound as the number of threads increases. The latter comes from an equilibrium between the threads added and the distribution of the workload. The distribution of work across threads tends to offset the increase in the number of thread pairs exactly and, consequently, the number of conflicts. 3.1.1. Validation on the DiSCo phantom To assess the effects of ExaTract parallelization, we run ExaTract on the DiSCo phantom using four parallelization strategies and compare with a baseline computed without parallelization. The effects of multithreading were evaluated on the non-partitioned domain using 16 threads in parallel. The effects of an extreme domain decomposition were isolated using a 16-partitioned domain (1 thread per partition) with ω = 10 op/voxel. Two intermediary experiments were conducted using 2 and 4 partitions, with 8 and 4 threads per partition, respectively. Global tractography is a stochastic optimization such that a single run differs from another when initialized with different seeds. To compare the variance due to the algorithm’s own randomness with the effects of parallelization, we computed 10 seed-independent runs per parallelization configuration. The obtained tractograms, compared qualitatively in Figure 6 and quantitatively in Table 2 in terms of fiber and spin-glass counts, show very similar fiber configurations, even for the extreme case of domain decomposition.

13

1 partition (1 thread)

1 partition (16 threads)

2 partitions 16 threads

4 partitions 16 threads

16 partitions 16 threads

Figure 6: Effects of the parallelization strategies on the DiSCo phantom using different domain decompositions and multithreading approaches while using a single seed per computation. Regardless of the decomposition scheme (top) and the computational resources allocated for multithreading, the resulting tractograms (middle) and their corresponding structural connectivity matrices (bottom) remain consistent.

In Table 2, we additionally report conflict occurrences. Conflict occurrences are consistent with previous estimates, with less than 0.07% of total operations in conflict in the extreme case of domain decomposition. As previously modeled, the multithreading conflicts are negligible compared to the communication conflicts, with a maximum of 9.73x10−5 % of the operations being in conflict for the extreme case of multithreading. We note that conflict occurrences vary from one experiment to another while using 16 threads in total (independently of the partitioning). It might come from the extreme cases of the parallelization on a very small dataset, where 1) overlapping areas are likely to be distributed over more than 2 MPI processes, which increases the probability of communication conflicts, and 2) the reduced dataset size leads to a higher probability for two threads to be synchronized (so that assumptions made in Appendix A are not valid) and draw to the same spin-glass, which increases the probability of multithreading conflicts. Despite these discrepancies, the observed conflicts remain negligible and do not significantly affect the system’s convergence. Since the DiSCo phantom is a synthetic dataset, the obtained tractograms were also compared to the ground truth using ROC AUC, which evaluates both tractogram sensitivity and accuracy (Girard et al., 2023). The ROC AUC stability indicates a strong reproducibility in the inferred tractograms across parallelization strategies. With a maximum standard deviation of 0.005, the ROC AUC distributions per experiment illustrate the robustness of the parallelization strategy to find consistent results from different initializations. To further compare the seed-induced variability to the parallelization-induced variability, we computed the Pearson correlations r seed to seed between the connectivity matrices obtained with different seeds in the same parallelization configuration and the Pearson correlations r parallel to sequential 14

between the connectivity matrices obtained in a given parallelization strategy to the baseline connectivity matrices computed sequentially. r seed to seed is very high, with a minimum value of 0.9977. This indicates that the variability induced by the stochastic nature of the algorithm is low and that the obtained tractograms are consistent across different runs. The comparison of the tractograms obtained with different parallelization strategies to the baseline tractogram computed without parallelization also shows a very high Pearson correlation r parallel to sequential , with a minimum value of 0.9978. This indicates that the parallelization strategies do not significantly affect the resulting tractograms and that the obtained results are consistent across different parallelization configurations. To assess whether partioning and multithreading parallelizations alter the connectivity beyond the inherent random floor of the global tractography optimization, we compared pairwise connectivity matrix correlations within the experiments of the same implementation (different seeds) and pairwise connectivity matrix correlations between the sequential and the parallel implementations. Because pairwise comparisons drawn from a shared pool of runs do not follow the standard assumptions of independent samples, we used a label-permutation test (10,000 resamples), which found no evidence that parallelization introduces divergence beyond the algorithm’s own noise: as reported in Table 2, all p-values were found greater than 0.05, rejecting the null hypothesis that the two distributions are different. Partitioning Threads per partition Overlap coverage (%) Fiber count Spin-glass count Multithreading conflicts (x10−7 ) Communication conflicts (x10−4 ) ROC AUC r seed to seed (45 pairs) r parallel to sequential (100 pairs) Permutation test p-value

1 1 0 26662 ±87 548875 ± 1223 0.0 0.0 0.0 ± 0.0 0.922 ± 0.003 0.9977 ±0.0005 -

1 16 0 26653 ±73 548470 ± 1450 0.39 ± 0.17 0.0 ± 0.0 0.923 ± 0.003 0.9980 ±0.0004 0.9978 ±0.0004 0.331

2 8 7 26789 ±78 550974 ± 1089 4.49 ± 0.67 3.79 ± 0.14 0.921 ± 0.005 0.9982 ±0.0003 0.9979 ±0.0004 0.408

4 4 18 26792 ±105 551556 ± 1582 9.73 ± 1.07 4.70 ± 0.09 0.922 ± 0.004 0.9979 ±0.0005 0.9978 ±0.0005 0.147

16 1 47 26906 ±77 551517 ± 899 0.0 ± 0.0 6.13 ± 0.05 0.922 ± 0.005 0.9980 ±0.0004 0.9979 ±0.0004 0.940

Table 2: Statistics of the ExaTract benchmark on the DiSCo dataset. To assess reproducibility, all experiments were repeated 10 times. Multithreading conflicts are reported as the ratio of corrupted operations out of the total operations performed. Communication conflicts are reported as the ratio of corrupted operations among total operations performed in overlapping regions.

3.1.2. Validation on an in vivo HCP human subject Similarly, five ExaTract runs were launched on one subject of the HCP database with 1, 8, 32, 64, and 128 partitions, multithreaded over 16 bi-threaded cores per partition. The obtained tractograms were automatically labeled, resampled in the MNI space, to extract the deep WM 15

bundles based on the deep WM bundle atlas established in Herlin et al. (2024). As compared in Figure 7, the shapes of the labeled bundles are similar from one run to another, even in the extreme case of 128 partitions. The obtained tractograms were also compared in terms of fiber and spin-glass count, and conflict occurrences (Table 3). The number of fibers and spin-glasses remains stable with the number of partitions, which is consistent. The occurrence of communication conflicts is consistent with the previous estimate, and is almost constant over the experiment, with a maximum of 0.023% of the operations in conflict for the extreme case of domain decomposition. Despite these conflicts, the obtained bundles are qualitatively similar, suggesting that the system’s convergence is not significantly affected. Regarding computation time, a significant acceleration is observed when increasing the number of partitions, with a speed-up of 88.4 using 128 partitions compared to the single-partition case. The parallelization efficiency remains high, at 0.69 for 128 partitions with overlapping regions covering 47% of the domain, indicating that the communication and synchronization overhead remains acceptable for the considered workload.

16

Domain decomposition

N=8

N = 128

Frontal aslant tracts

IFOF

Cortico-spinal tract

Arcuate bundle

N=1

Figure 7: Fiber labeling results on one subject of the HCP database from three ExaTract runs computed with 32 (N = 1 partition), 256 (N = 8 partitions), and 4096 (N = 128 partitions) threads, on the arcuate bundle, the corticospinal tracts, the inferior fronto-occipital fasciculus, and the frontal aslant tracts. The shape of the recovered bundles is similar from one run to another.

17

Partitioning Overlap coverage (%) Fiber count Spin-glass count Multithreading conflicts (x10−7 ) Communication conflicts (x10−4 ) Computation time Speed-up Efficiency

1 0 886186 2.94x107 2.35x10−3 0 17h24 1 1

8 9.2 889966 2.96x107 1.78 1.79 2h40 6.5 0.82

32 21 891224 2.97x107 4.51 2.08 44min35s 23.4 0.73

64 31.88 894942 2.98x107 6.48 2.09 22min54s 44.6 0.71

128 45.1 899740 3.0x107 9.11 2.27 11min48s 88.4 0.69

Table 3: Tractograms statistics, conflict occurrences, and performance of ExaTract benchmarked on one subject of the HCP dataset. 16 bi-threaded cores are used for each partition.

3.2. Scalability benchmarking at ultra high resolution The scalability of ExaTract was evaluated under two perspectives: strong and weak scaling. In the first case, scalability is evaluated by measuring parallelization efficiency on a fixed-size problem, while in the second, it is evaluated on a growing dataset while preserving the workload per processor. To evaluate the strong scaling, ExaTract was run on the ex vivo macaque brain dMRI dataset using increasing numbers of cores while keeping the entire tractography domain fixed. For each combination, the speed-up of the first 20 ExaTract iterations was measured for three different workloads (ω = 2/5/20 op/voxel) and normalized to baselines. The respective baselines were computed from the same workloads performed using 16 cores only. The macaque dMRI dataset was chosen for this experiment because it is the largest dataset that can be processed on a single computing node of the Irène-Rome cluster, allowing us to evaluate strong scaling without introducing additional overhead from inter-node communication. The strong scaling demonstrates the benefits of parallelization, achieving reduced computation time by effectively reducing the workload per processor without introducing excessive overhead. The weak scaling curve was established using the post-mortem human brain Chenonceau dataset to illustrate scaling behavior on an artificially growing dataset while preserving the workload per processor. To this aim, the tractography domain of the Chenonceau brain was split into 800 parcels. From 128 cores to 12800 cores, corresponding to 1 and 100 nodes, respectively, 8 to 800 parcels were involved, using 16 cores per MPI process. Similar to the strong scaling curve, the speed-up was normalized to the time required for the same workload performed on a single parcel. As expected, the workload plays an important role in the performance. While a heavy workload improves computational performance, a low workload requires more communication to reduce errors. Both experiments show correct scalability (Figure 8). While the strong scaling curve shows significant acceleration on a given dataset, the weak scaling curve demonstrates ExaTract’s ability to handle very high-resolution, large datasets that require substantial memory. The linear asymptotic behavior of the weak scaling curve suggests that scalability is even feasible for larger datasets. This property is strategically aligned with the emergence of very large datasets in microscopy, combining dMRI, 3D-PLI, or comSLI, such as the BigMac dataset (Howard et al., 2023). Thanks to the assumptions made, especially the domain decomposition, an MPI process dealing with a given brain parcel interacts only with the processes in charge of neighboring parcels, ignoring the rest. It means upscaling depends only on the available resources for the largest tractography domains, which supports high scalability potential on future 18

exascale supercomputers. Strong scaling

Weak scaling

Ideal speed-up =2 =5 = 20

Ideal speed-up =1 =5

104

Speed-up

104

103 103

81 92 12 80 0

Cores number

40 96

20 48

10 24

51 2

25 6

12 8

8

4

32 76

Cores number

16 38

81 92

40 96

20 48

10 24

51 2

25 6

102

Figure 8: Scalability curves of ExaTract on the Irène-Rome CPU partition. The speed-up was evaluated in strong (left) and weak (right) scaling problems using from 1 node (128 cores) up to 256 nodes of the Irène-Rome partition.

3.3. Proof of scalability for population imaging Before testing ExaTract against the entire HCP dataset, we benchmarked the performance of the designed parallelized tractography algorithm on the available hardware by mapping the performance landscape across various deployments on the Irène-Rome computer. The tractography domain of one HCP subject was partitioned into powers of 2 from 1 to 256. For each partitioning, the number of threads per MPI process and the required number of nodes (up to one partition per node) were varied. Different tractography experiments were run with an evaluation of the first 10 MPI iterations’ performance. The speed-up and corresponding efficiency are respectively shown in Figure 9. The optimal settings that yield the best performance at fixed numbers of cores are represented by the blue curve in Figure 9. Up to 2048 cores (corresponding to 16 nodes), the performance profile, corresponding to these optimal settings (blue curve), is also reported in Figure 10.

19

Acceleration

Efficiency 1.0 1000

0.8

800

1.0

600

0.6

400

0.2

0.8

1200 1000 800 600 400 200

0.4

sed

128

64

32 Partio 16 8 nning

12

200

cor es u

256

32768 16384 8192 4096 2048 1024 512 256 128 64 32 16 8 4 2 1

4

2

1

48 1632 cor 64128 es 256 use 512 d 1024 2048 4096 8192 16384 32768

64

1

2

4

8

16

128

256

0.6

0.4

0.2

32 ing

onn

ti Par

Efficiency

Optimal threads count per MPI rank 30 25 Thread count

20 15 10 5

Number of CPU cores

1 2 4 8 16 32 64 12 8 25 6 51 10 2 2 20 4 48

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3

1 2 4 8 16 32 64 12 8 25 6 51 10 2 2 20 4 48

Efficiency

Figure 9: Left: Speed-up landscape across different deployment configurations on the Irène-Rome CPU partition. The speed-up is calculated using the reference time obtained with one MPI process launched on one core. The blue curve illustrates the optimal configuration that provides the fastest tractography at a fixed number of cores. Right: efficiency across different deployment configurations. The reference time stemming from a deployment on one MPI process and 1 core was used to compute the efficiency.

Number of CPU cores

Figure 10: Left: profile of the optimal efficiency obtained among the tested deployments on an HCP subject. Right: corresponding optimal number of threads used per MPI process (or per partition). Over 16 cores used, the optimal number of threads used per MPI process never exceeds 32, illustrating the need to split the global tractography task into different MPI processes inside a node.

At millimeter resolutions, the scalability up to 32,768 cores (or 256 nodes) is, as expected, nonsensical. The experiment was primarily conducted for ExaTract behavior characterization. The efficiency remains above 50% up to 256 cores (2 nodes), whereafter it drops below this limit. Interestingly, the optimal configurations, as displayed in Figure 9 and Figure 10, show that the optimal number of threads per MPI process is never the absolute maximum of 256 but 20

varies between 8 and 32. This is due to the non-uniform memory accesses of the NUMA nodes of the Irène-Rome computers. Memory access times, therefore, depend on the thread’s position and the location of the allocated memory. It plays in favor of the use of more MPI processes to better benefit from the use of memory proximity tailored to the hardware. For this reason, the number of threads per MPI process was fixed at 32, corresponding to 16 cores, for the following experiments on the HCP and Chenonceau datasets. Due to the Irène-Rome supercomputer scheduler’s preference for whole-node usage, we constrained the ExaTract deployment to a single HCP subject per NUMA node. In that case, the optimal configuration for 128 cores is to use 8 MPI processes, each with 32 threads. Using 178 nodes of the Irène-Rome HPC facility simultaneously (totaling 22784 CPU cores), the application of the ExaTract tractography method over the 1065 subjects took less than 24 hours. With 8.07x105 streamlines generated on average per subject in 3 h 40 min, the obtained tractograms are at least as comprehensive in terms of fiber count as those generated by the conventional streamline algorithm using one seed per voxel to create a streamline. The correspondence between the computation time per subject and its tractography domain size is shown in Figure 12. In total, 8.59x108 fibers were generated across the 1065 subjects in 22 hours. This makes it the largest global tractography dataset ever computed, in terms of fiber count. A small sample of these results is displayed in Figure 11.

21

Figure 11: Coronal views of 20 randomly selected tractograms in their native space computed with ExaTract among the 1065 subjects of the HCP cohort.

22

1e6 1.4

6 5

1.2 1.1

Generated fiber count

ExaTract computation time (h)

1.3

4

1.0

4

0.9 4

0.8 0.7

3 500,000

600,000

700,000

Voxel count

800,000

900,000

0.6

Figure 12: Global tractography Computation time with ExaTract for each HCP subject with respect to voxel count in the tractography domain, colored by the number of generated streamlines. The plot shows the expected linear relationship between the two variables. The voxel count directly determines the number of operations, and with a similar density of spin-glasses per voxel, the computation time increases linearly with the tractography domain size.

3.4. Proof of scalability on an ultra-high resolution dataset ExaTract is not only a valuable tool for population imaging analysis; it is also of major interest in the community to deal with ultra-high mesoscopic resolution datasets. In this second application, we focused on the use of ExaTract to infer the structural connectivity of the Chenonceau dataset, an ultra-high-field diffusion MRI dataset of a post-mortem human brain scanned at the mesoscale. The tractography domain of the Chenonceau brain was partitioned into 800 parcels (see Figure 13A) that were distributed over 100 nodes of the Irène-Rome CPU partition, corresponding to 12800 bi-threaded CPU cores. Given the previous benchmark of section 3.2, the efficiency of this deployment is expected to be around 0.6 at ω = 5 op/voxel compared to a 32-threaded approach, which remains a good trade-off between workload and communication. The number of threads per MPI process was set to 32 to favor local memory usage. Three runs were performed to reach convergence, with a gradual increase in the workload between runs and a decrease in the cooling schedule to reach convergence. At the end of each run, temperature schedules were empirically adapted to keep the convergence tendency. The three simulated annealing cycles comprised 10,000 MPI iterations each with an exponential decrease from 0.2 to 0.003 K for the first cycle, from 0.1 to 0.003 K for the second cycle, and from 0.05 to 0.003 K for the third cycle (see Figure 13B). Initial spin-glass per voxel was set to ninit = 4 spin-glass/voxel and l = r search = 160 µm. To speed up convergence, the number of sampled directions in the ODFs was reduced to 1000, and ω was increased from 5 to 10 op/voxel for the second and third runs. In Figure 13C-F, the timing and convergence analysis of the three runs are displayed. The optimization kernels, which represent the core of the tractography process, are predominant compared to the expansion and consistency kernels. Communication overhead remains acceptable compared to the optimization time and in agreement with the previous benchmark. It is worth noting that the regular peaks observed in the communication kernel correspond 23

to intermediate saves of the spin-glass map to the storage disk (every 500 iterations). The ratio of accepted random motion is highly dependent on the system temperature, allowing for reconfigurations and exploration by altering spin-glass positions and orientations. The evolution of the accepted modification ratio for creation, deletion, and connection operations illustrates the system’s expected convergence, with a continuous decrease in acceptance rates throughout the progression of the tractography. Similarly, the number of spin-glasses and free extremities across the three runs correlates with the expected behavior, with a monotone decrease of free extremities and a logarithmic growth of the number of spin-glasses reaching an asymptote, reflecting the direct effects of the connection, deletion, and creation operations. Conflict occurrence remains very low relative to the total number of operations and spin-glasses, validating the proposed parallelization scheme.

24

B

A

Cooling schedule 0.200 0.175

temperature (K)

0.150 0.125 0.100 0.075 0.050 0.025 0.000 0

C

5000

D

Kernels timing

10000

15000 iterations

20000

25000

30000

25000

30000

Spin-glass and extremity count 107

106

count

time (ms)

106

105

105 communication optimization expansion consistency

104 0

5000

E

10000

15000 iterations

20000

25000

30000

0

5000

F

Conflict occurence

10000

15000 iterations

20000

Accepted operations 100

domain decomposition multithreading

10

acceptance ratio

102

conflicts count

spin-glass free extremity

104

101

motion connection creation deletion

1

10 2

10 3

10 4 100 0

5000

10000

15000 iterations

20000

25000

30000

0

5000

10000

15000 iterations

20000

25000

30000

Figure 13: Application of ExaTract to the Chenonceau dataset. (A) The entire brain was partitioned into 800 parcels using METIS (Karypis and Kumar, 1998). (B) The cooling schedule was tuned to ensure convergence over three runs, with a logarithmic temperature decrease. (C-F) Timing and convergence analysis performed on rank 0. (C) Kernel timings across ExaTract iterations reveal the dominance of optimization kernels over expansion and consistency kernels. Communication overhead remains comparable to the optimization cost, which is an expected trade-off when scaling to datasets of this size. (D) The evolution of spin-glass count and free extremities across the three runs directly confirms the system’s proper convergence. (E) Conflict occurrence remains negligible relative to the total number of operations. (F) The evolution of the accepted modification ratio across operation types further illustrates the system’s reconfiguration dynamics during the convergence.

It respectively took 18.6, 16.3, and 17.9 hours for the three runs to be completed, totaling 52.8 hours of consecutive computation. At a machine level, it represents 675,840 core hours of cumulative computation. A total of 4.2x109 spin-glasses were involved in the 1.74x1013 distributed operations, generating 25x106 fibers. Figure 14 illustrates the obtained tractogram computed from ExaTract, with an overview of anatomical details obtained thanks to the combination of the global tractography algorithm with the high-spatial and angular resolutions of the Chenonceau dataset. To further assess the correct convergence of the tractography, the deep WM bundles were automatically labeled similarly to the methods developed in section 3.1.2. The obtained bundles are displayed in Figure 15 and show a good recovery of the well-established deep WM bundles, such as the arcuate fasciculus, the inferior fronto-occipital fasciculus, the inferior lon25

gitudinal fasciculus, the uncinate fasciculus, the cingulum bundle, the corpus callosum, and the corticospinal tract. As presented in Legeay et al. (2026), the tractogram was further processed with clustering methods to extract fiber bundles. Figure 15 depicts all bundles according to their length range. The amount of found bundles using global tractography exposes the richness of the short connectivity in the human brain, from very short intra-cortical connections up to deeper bundles connecting adjacent gyri.

26

Figure 14: Details of the global tractography results obtained on the Chenonceau brain in its native space. The entire tractogram yields approximately 25 million fibers; here, 1/4 of the total fibers are displayed in the axial view (A). Fibers in a crossing area (B) show the ability of the global tractography to recover complex fiber configurations, such as crossing. When zooming on a sulcus (C), superficial WM bundles joining two gyri are immediately visible (white dotted lines). Deep in the sulcus, fibers in the gray matter keep a radial trajectory (blue lines). Interestingly, some fibers originating in the cortex run along the gyrus (orange arrows), indicating even shorter connections.

27

Figure 15: Global tractography atlas of the Chenonceau brain. Fiber bundles stemming from a clustering of the tractogram (see Legeay et al. (2026)) were sorted per length to28 expose the different connectivity levels. Very short connections, ranging here from 5 to 25 mm, remain particularly superficial, or even intra-cortical, while longer connections extend deeper into the white matter. The established deep WM bundles were extracted using an automatic procedure based on the atlas of Herlin et al. (2024). It includes the inferior longitudinal fasciculi, the thalamic radiations, split into the temporal (turquoise blue), frontal (blue), central (orange), parietal (pink) and occipital (green) thalamic radiations, the frontal aslant tracts, the frontal (green), central (blue), parietal (turquoise blue), pink (red) lenticular radiations, the cortico-spinal tracts (CST), the corpus callosum, the inferior fronto-occipital fasciculi (IFOF), the fornix, the cingulum bundles, split into the dorsal, long (red) and short (blue), and dorso-ventral (green), the middle longitudinal fasciculi, the uncinate bundles, the SLF, including the SLF I (blue), II (yellow) and III (purple), the optic radiations, the cingulum ventral, the arcuate bundles.

4. Discussion Thanks to spatial decomposition and an efficient multithreading strategy, ExaTract effectively reduces computation time while allowing more fibers to be generated. Through dedicated experiments on the DiSCo phantom and an HCP subject, we demonstrated that the parallelization scheme does not significantly affect the system’s convergence, even in extreme domaindecomposition cases. The obtained tractograms remain similar across different parallelization strategies, and the conflict occurrences remain negligible compared to the total number of operations performed. It remains, however, challenging to further characterize convergence in overlapping areas, as conflicts may have additional effects on the system’s convergence. While the experiments do not show significant differences in the resulting tractograms, it is possible that convergence is locally and dynamically affected in overlapping areas, which may not be visible at the global level or at the end of the convergence process. Future work could focus on further characterizing the convergence in these areas dynamically, yet, due to the randomness of the MCMC process and the many possible variations in the energy landscape, these effects might be difficult to characterize and may require a very large number of experiments to be statistically significant. The scalability benchmarks further illustrate ExaTract’s ability to effectively reduce computation time while maintaining high efficiency, even at large numbers of cores. In comparison, tools like MRtrix (Reisert et al., 2011; Christiaens et al., 2015) or MITK (Neher et al., 2012) face time and memory constraints that limit the number of generated fibers (typically < 100,000) and the resolution of the datasets they can handle. This can fall short when representing the structural connectivity of the entire human brain. We reported in Table 4 the performance of different global tractography tools on human datasets. While the comparison is difficult due to differences in the experimental setup, it provides an order-of-magnitude estimate of ExaTract’s performance relative to existing tools. In a similar computational time, ExaTract generates 10 times more fibers than the tool from the previous implementation, while optimizing 10 more spin-glasses in a slightly more resolved dataset. As shown earlier, with 4096 threads, ExaTract achieved the same fiber-generation performance as the 256-thread implementation in less than 12 minutes, making it suitable for clinical research investigations. It is worth noting that the comparisons remain highly sensitive to the experimental setup: first, computing hardware design technologies have advanced substantially since the first demonstration of global tractography. Second, the chosen initial conditions of the global tractography algorithm can significantly impact the global tractography process. The initial spin-glass density per voxel directly influences the total number of generated streamlines and increases the workload when searching for potential connections. Other parameters linked to the diffusion model, such as the number of directions used to sample the ODF or the choice of a generative spin-glass model, directly affect the workload. Finally, spin-glass length and research diameter can affect the convergence speed. In addition to time constraints, current implementations cannot handle very large datasets, such as the Chenonceau dataset, due to memory limitations. ExaTract, on the other hand, is designed to handle such large datasets by distributing the workload across multiple nodes and effectively managing memory usage through domain decomposition. This allows ExaTract to perform global tractography on ultra-high-resolution datasets that would be infeasible with stateof-the-art tools. An additional advantage of ExaTract is its available memory obtainted through spatial decomposition. Computation at the level of the diffusion model used is often demanding due to the limited memory available in previous implementations of tractography. With addi29

tional memory, ExaTract could keep the orientation probabilities in RAM for each voxel, which are currently computed on the fly. This would allow us to save a lot of time by avoiding the repeated computation of the diffusion model for each spin-glass and each operation. Global tractography implementation Fillard et al. (2009) Reisert et al. (2009) Wu et al. (2016) Konopleva et al. (2018) ExaTract (HCP) ExaTract (HCP) ExaTract (Chenonceau)

resolution (mm3 ) 1.8x1.8x2 2x2x2 2x2x2 1.5x1.5x3 1.25x1.25x1.25 1.25x1.25x1.25 0.2x0.2x0.2

spin-glass count 2x106 2.5x106 2.5x106 3x106 3x107 3x107 4.2x109

fiber count 105 2x105 8x105 8x105 2.5x107

computation time 3 days <24h 1h30min 6h 2h40min 11min48s 52h48min

threads used 10 2 256 4096 12800

Table 4: Comparison of fiber count and convergence timings for multiple tractography implementations in human datasets.

A strong limitation of global tractography is the difficulty in predicting the convergence time. The number of generated fibers and spin-glasses, as well as the convergence speed, are highly dependent on the initial conditions of the tractography. This makes it challenging to estimate the computational resources and time required for a given dataset, especially for large, highresolution datasets. The additional parameters involved in the parallelization add complexity to setting up the tractography. Users need to find the right balance between workload and communication, which can be difficult without prior experience or benchmarks. This can lead to suboptimal performance if not properly configured. Future work could focus on developing more adaptive parallelization strategies that dynamically adjust the workload and communication based on dataset characteristics and available computational resources, potentially improving efficiency and reducing the time required for global tractography. One of the key advantages of the portability library Kokkos used in ExaTract is its capability to support emerging hardware architectures with minimal code modifications, from CPU to GPU-based architectures. While a non-generative approach for the spin-glass model was implemented in this work, ExaTract is designed to be flexible and adaptable to various energy models, opening the way to plug any local model of the diffusion process and any regularization term to improve the robustness of the tractography. This flexibility allows ExaTract to be used with different global tractography frameworks. However, generative approaches, as first introduced by Kreher et al. (2008), involve a modeling of the contribution of each spin-glass in the diffusion-weighted signal. The optimization process seeks to minimize the voxel-wise difference between the generated diffusion signal from spin-glasses and the actual diffusion data by expressing the external energy as follows: 2

Eext (M , D) =

X

X

ρ Xi − D V

,

(8)

voxel V Xi ∈V

where ρXi is the signal generated by the spin-glass Xi according to a modelling assumption and DV is the measured diffusion data in voxel V. Reisert et al. (2014), Christiaens et al. (2015) and Konopleva et al. (2018) have recently unified diffusion modeling and global tractography. Under the name of model-free global tractography, Konopleva et al. (2018) makes the minimal assumption possible on the fiber response function of a spin-glass, meaning that the fiber response is 30

re-estimated per voxel after each modification of the spin-glass configuration by an optimization under constraint over the fiber response function. In these cases, the external energy is computed voxel-wise, with all spin-glasses in each voxel taken into account to estimate the model. This approach requires some adaptation of ExaTract’s parallelization scheme, since conflict occurrences are likely to increase when performing the reduction at the voxel level. A possible adaptation of the model-free global tractography to ExaTract could be to perform the reduction at the iteration level, meaning that the fiber response function is updated at the end of each MPI iteration, thereby allowing the parallelization scheme to remain unchanged. 5. Conclusion We demonstrate the scalability and performance of ExaTract on large datasets, showcasing its ability to reconstruct brain connectivity at unprecedented scales while maintaining accuracy. The entire HCP dataset was processed in under 24 hours using ExaTract, demonstrating its capability to handle large-scale neuroimaging data. The application of ExaTract to the ultra-high-resolution Chenonceau dataset further illustrates its ability to reconstruct detailed brain connectivity patterns at a mesoscopic scale. This achievement allowed us to build the first mesoscopic structural connectome from global tractography. By addressing the computational challenges associated with global tractography, ExaTract opens new avenues for exploring the advantages of global approaches in both research and clinical settings. ExaTract is completely open-source and available on Framagit as part of the ginkgo toolkit. Acknowledgments This research was supported by the Joint German-French Institute AIDAS Cycle 1 and Cycle 2, by the European Union’s Horizon 2020 Framework Program for Research and Innovation under the Specific Grant Agreement No. 945539 (Human Brain Project SGA2 and SGA3), and by the Priority Research and Equipment Program (PEPR) Santé Numérique BrainDeepPhenotyping. All the acquisitions performed for the Chenonceau dataset were entirely covered by the NeuroSpin imaging platform. This project was provided with computing HPC and storage resources by GENCI at TGCC, thanks to the grants A0120310851 and SS010315366 on the supercomputer Joliot Curie’s ROME partition. The authors would also like to thank Dr. France Boillod-Cerneux and Dr. Christophe Calvin for their support within the framework of the French-German AIDAS institute, as well as the GENCI French national HPC infrastructure. Data were provided in part by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657), funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. Appendix A. Conflict occurence estimation In the case of parallelization on a shared-memory architecture, the probability of conflict occurrence can be estimated as follows. During an ExaTract iteration, the total workload performed on a partition with Nv voxels, with n spin-glasses each, is divided by all Nt working threads, so that each thread performs NNv ωt operations per iteration. For a given spin-glass, the probability to 31

Nv ω

be drawn by a thread is 1 − (1 − N1 ) Nt . Supposing that all operations have a duration t and that the beginning of operations on two threads are independent and uniformly distributed in [0, t], the probability of overlapping modifications for two threads is 1 − (1 − ttc )2 , where tc is the critical zone duration. Hence, the probability of multithreading conflict occurrence X for a given spin-glass can be estimated as follows:  ! Nv ω 2 !   tc 2 2 1 Nt    1 − 1 − R p(X) = 1 − 1 − nNv t

(A.1)

  Given the N2t different pairs of threads that can potentially modify the same spin-glass, and according to the linearity of the expectation over the entire spin-glass set, the mean conflict occurrence during an iteration, normalized per operation, is given by: ! E[X] Nt n p(X) (A.2) = ωNv 2 ω Supposing that modification draws are independent and uniformly distributed per voxel, and that a modification made in the overlapping area of a partition has a probability of poverlap to alter a spin-glass in a neighboring partition, the expectation of communication conflict occurrence Y in an overlapping area can be estimated as follows: E[Y] 1 n 1− 1− = ωNv ω n

!ω !2

2

R poverlap

(A.3)

Appendix B. Conflict resolution strategy The approach for resolving conflicts is based on a kernel that detects inconsistencies in the spin-glass configurations and corrects them according to a predefined resolution plan. For this strategy to be effective, the merging of the modifications performed in the overlapping area was thought to be symmetrical and deterministic, ensuring that the same correction is consistently applied on all sides of the overlap without additional communication. This is achieved through a strategy of merging with a posteriori corrections, which is divided into the following four steps: 1. Identification of modifications: All modifications made within the overlap area are identified by comparing the current state with a saved copy of the overlap area before any modifications were performed. 2. Exchange of modifications: Using the communication protocol, all these modifications are exchanged among neighboring nodes that share the same overlap area. 3. Application of modifications: In a predefined order, arbitrarily chosen and planned in advance by all the nodes, the changes made by all the nodes sharing the overlap area are sequentially applied to the original copy of the overlap area. This process can lead to conflicts, as modifications could have been performed in parallel. Despite potential conflicts, all modifications are applied, potentially resulting in inconsistent spin-glass configurations, as illustrated in Figure B.16. 4. Detection and correction of abnormalities: Any abnormality, such as spin-glass inconsistencies within the overlap area, is detected and rectified according to a predefined resolution plan in which all possible configurations are anticipated. Figure B.16 exposes an example of a possible conflict resolution. 32

This approach enables rapid recovery from conflicts, efficiently restoring the system to a consistent state with a linear complexity kernel. Node A

Node B

Initial configuration

Parallel computation during an MPI iteration

Applying the modification from node A

Applying the modification from node B

Correction of the conflict according to a planned decision

Figure B.16: Example of conflict and subsequent correction. Node A and node B start from the same initial configuration in a voxel shared in the overlapping area between the two nodes. While node A connects the red spin-glass to the purple one, node B independently connects the red spin-glass with the orange one. After communicating their respective modifications, the modification from node A is first applied to the previously stored initial configuration. The modification from node B is then applied, leaving the purple spin-glass non-reciprocally connected to the red spin-glass. This leads to an anomaly that is further detected and corrected according to a predefined resolution plan.

References Axer, M., Grässel, D., Kleiner, M., Dammers, J., Dickscheid, T., Reckfort, J., Hütz, T., Eiben, B., Pietrzyk, U., Zilles, K., Amunts, K., 2011. High-Resolution Fiber Tract Reconstruction in the Human Brain by Means of Three-Dimensional Polarized Light Imaging. Front Neuroinform 5, 34. https://doi.org/10.3389/fninf.2011.00034 Axer, M., Amunts, K., 2022. Scale matters: The nested human connectome. Science 378, 500–504. https://doi.org/10.1126/science.abq2599 33

Carter Edwards, H., Trott, C.R., Sunderland, D., 2014. Kokkos: Enabling manycore performance portability through polymorphic memory access patterns. Journal of Parallel and Distributed Computing, Domain-Specific Languages and High-Level Frameworks for High-Performance Computing 74, 3202–3216. https://doi.org/10.1016/j.jpdc.2014.07.003 Christiaens, D., Reisert, M., Dhollander, T., Sunaert, S., Suetens, P., Maes, F., 2015. Global tractography of multi-shell diffusion-weighted imaging data using a multi-tissue model. NeuroImage 123, 89–101. https://doi.org/10.1016/j.neuroimage.2015.08.008 L. Clarke, I. Glendinning, and R. Hempel, The MPI Message Passing Interface Standard, in: K. M. Decker and R. M. Rehmann (Eds.), Programming Environments for Massively Parallel Distributed Systems, Monte Verità, Birkhäuser, Basel, pp. 213–218, 1994. Cointepas, Y., Mangin, J.-F., Poupon, C., 2002. A Spin Glass Based Framework to Reconstruct Brain Fiber Bundles from Images of the Water Diffusion Process. Daducci, A., Dal Palù, A., Lemkaddem, A., Thiran, J.-P., 2015. COMMIT: Convex Optimization Modeling for Microstructure Informed Tractography. IEEE Transactions on Medical Imaging 34, 246–257. https://doi.org/10.1109/TMI.2014.2352414 Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R., 2007. Regularized, fast, and robust analytical Q-ball imaging. Magnetic Resonance in Medicine 58, 497–510. https://doi.org/10.1002/mrm.21277 Fillard, P., Poupon, C., Mangin, J.-F., 2009. A novel global tractography algorithm based on an adaptive spin glass model. Med Image Comput Comput Assist Interv 12, 927-934. Girard, G., Rafael-Patiño, J., Truffet, R., Aydogan, D.B., Adluru, N., Nair, V.A., Prabhakaran, V., Bendlin, B.B., Alexander, A.L., Bosticardo, S., Gabusi, I., Ocampo-Pineda, M., Battocchio, M., Piskorova, Z., Bontempi, P., Schiavi, S., Daducci, A., Stafiej, A., Ciupek, D., Bogusz, F., Pieciak, T., Frigo, M., Sedlar, S., Deslauriers-Gauthier, S., Kojčić, I., Zucchelli, M., Laghrissi, H., Ji, Y., Deriche, R., Schilling, K.G., Landman, B.A., Cacciola, A., Basile, G.A., Bertino, S., Newlin, N., Kanakaraj, P., Rheault, F., Filipiak, P., Shepherd, T.M., Lin, Y.-C., Placantonakis, D.G., Boada, F.E., Baete, S.H., Hernández-Gutiérrez, E., RamírezManzanares, A., Coronado-Leija, R., Stack-Sánchez, P., Concha, L., Descoteaux, M., Mansour L., S., Seguin, C., Zalesky, A., Marshall, K., Canales-Rodríguez, E.J., Wu, Y., Ahmad, S., Yap, P.-T., Théberge, A., Gagnon, F., Massi, F., Fischi-Gomez, E., Gardier, R., Haro, J.L.V., Pizzolato, M., Caruyer, E., Thiran, J.-P., 2023. Tractography passes the test: Results from the diffusion-simulated connectivity (disco) challenge. NeuroImage 277, 120231. https://doi.org/10.1016/j.neuroimage.2023.120231 Green, P.J., 1995. Reversible jump Markov chain monte carlo computation and Bayesian model determination. Biometrika 82, 711–732. https://doi.org/10.1093/biomet/82.4.711 Greening, D.R., 1990. Parallel simulated annealing techniques. Physica D: Nonlinear Phenomena 42, 293–306. https://doi.org/10.1016/0167-2789(90)90084-3 Herlin, B., Uszynski, I., Chauvel, M., Dupont, S., Poupon, C., 2024. Sex-related variability of white matter tracts in the whole HCP cohort. Brain Struct Funct 229, 1713–1735. https://doi.org/10.1007/s00429-024-02833-0 34

Howard, A.F.D., Huszar, I.N., Smart, A., Cottaar, M., Daubney, G., Hanayik, T., Khrapitchev, A.A., Mars, R.B., Mollink, J., Scott, C., Sibson, N.R., Sallet, J., Jbabdi, S., Miller, K.L., 2023. An open resource combining multi-contrast MRI and microscopy in the macaque brain. Nat Commun 14, 4320. https://doi.org/10.1038/s41467-023-39916-1 Jayaraman, R., Darema, F., 1988. Error tolerance in parallel simulated annealing techniques. Presented at the Proceedings 1988 IEEE International Conference on Computer Design: VLSI, IEEE Computer Society, p. 545,546,547,548-545,546,547,548. https://doi.org/10.1109/ICCD.1988.25759 Karypis, G., Kumar, V., 1997. METIS: A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices (Report). Karypis, G., Kumar, V., 1998. Multilevel Algorithms for Multi-Constraint Graph Partitioning, in: Proceedings of the IEEE/ACM SC98 Conference. Presented at the SC98 - High Performance Networking and Computing Conference, IEEE, Orlando, FL, USA, pp. 28–28. https://doi.org/10.1109/SC.1998.10018 Kjolstad, F.B., Snir, M., 2010. Ghost Cell Pattern, in: Proceedings of the 2010 Workshop on Parallel Programming Patterns. Presented at the ParaPLoP ’10: 2010 Workshop on Parallel Programming Patterns, ACM, Carefree Arizona USA, pp. 1–9. https://doi.org/10.1145/1953611.1953615 Konopleva, L., Il’yasov, K.A., Skibbe, H., Kiselev, V.G., Kellner, E., Dhital, B., Reisert, M., 2018. Modelfree global tractography. NeuroImage 174, 576–586. https://doi.org/10.1016/j.neuroimage.2018.03.058 Kreher, B. w., Mader, I., Kiselev, V. g., 2008. Gibbs tracking: A novel approach for the reconstruction of neuronal pathways. Magnetic Resonance in Medicine 60, 953–963. https://doi.org/10.1002/mrm.21749 Kruper, J., Rokem, A., 2026. Massively parallelized brain tractography using compute clusters, supercomputers, and graphics processing units. Imaging Neuroscience. https://doi.org/10.1162/IMAG.a.1341 Legeay, S., Beaujoin, J., Hot, R.Y., Popov, A., Uszynski, I., Herlin, B., Poupon, F., Maldonado, I.L., Destrieux, C., Poupon, C., 2026. Chenonceau: mesoscopic MRI deep phenotyping of a post-mortem human brain at 7 and 11.7 Tesla. BioRXiv. https://doi.org/10.64898/2026.05.26.726838 Mangin, J.-F., Poupon, C., Cointepas, Y., Rivière, D., Papadopoulos-Orfanos, D., Clark, C.A., Régis, J., Le Bihan, D., 2002. A framework based on spin glass models for the inference of anatomical connectivity from diffusion-weighted MR data – a technical review. NMR in Biomedicine 15, 481–492. https://doi.org/10.1002/nbm.780 Mangin, J.-F., Fillard, P., Cointepas, Y., Le Bihan, D., Frouin, V., Poupon, C., 2013. Toward global tractography. NeuroImage, Mapping the Connectome 80, 290–296. https://doi.org/10.1016/j.neuroimage.2013.04.009 35

Menzel, M., Reuter, J.A., Gräßel, D., Huwer, M., Schlömer, P., Amunts, K., Axer, M., 2021. Scattered Light Imaging: Resolving the substructure of nerve fiber crossings in whole brain sections with micrometer resolution. Neuroimage 233, 117952. https://doi.org/10.1016/j.neuroimage.2021.117952 Neal, R.M., 1993. Probabilistic Inference Using Markov Chain Monte Carlo Methods. Neher, P., Stieltjes, B., Reisert, M., Reicht, I., Meinzer, H.-P., Maier-Hein, K., 2012. MITK global tractography. Presented at the Proceedings of SPIE - The International Society for Optical Engineering, p. 83144D. https://doi.org/10.1117/12.911215 Neiswanger, W., Wang, C., Xing, E.P., 2014. Asymptotically exact, embarrassingly parallel MCMC, in: Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, UAI’14. AUAI Press, Arlington, Virginia, USA, pp. 623–632. Philips, S., Hlawitschka, M., Scheuermann, G., 2013. Parallelized Global Brain Tractography. Vision. https://doi.org/10.2312/PE.VMV.VMV13.097-104 Poupon, C., Clark, C.A., Frouin, V., Régis, J., Bloch, I., Le Bihan, D., Mangin, J.-F., 2000. Regularization of Diffusion-Based Direction Maps for the Tracking of Brain White Matter Fascicles. NeuroImage 12, 184–195. https://doi.org/10.1006/nimg.2000.0607 Poupon, C., Destrieux, C., Maldonado, I. L., Uszynski, I. (2025). High-resolution postmortem anatomical ultra-high field MRI of the macaque brain (v1) [Data set]. EBRAINS. https://doi.org/10.25493/S5AC-HGG Rafael-Patino, J., Girard, G., Truffet, R., Pizzolato, M., Caruyer, E., Thiran, J.-P., 2021. The diffusion-simulated connectivity (DiSCo) dataset. Data in Brief 38, 107429. https://doi.org/10.1016/j.dib.2021.107429 Reisert, M., Mader, I., Kiselev, V., 2009. Global reconstruction of neuronal fibres. MICCAI Workshop on Diffusion Modelling. Reisert, M., Mader, I., Anastasopoulos, C., Weigel, M., Schnell, S., Kiselev, V., 2011. Global fiber reconstruction becomes practical. NeuroImage 54, 955–962. https://doi.org/10.1016/j.neuroimage.2010.09.016 Reisert, M., Kiselev, V.G., Dihtal, B., Kellner, E., Novikov, D.S., 2014. MesoFT: Unifying Diffusion Modelling and Fiber Tracking, in: Golland, P., Hata, N., Barillot, C., Hornegger, J., Howe, R. (Eds.), Medical Image Computing and Computer-Assisted Intervention - MICCAI 2014. Springer International Publishing, Cham, pp. 201-208. Smith, R.E., Tournier, J.-D., Calamante, F., Connelly, A., 2015. SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography. NeuroImage 119, 338–351. https://doi.org/10.1016/j.neuroimage.2015.06.092 Sudlow, C., Gallacher, J., Allen, N., Beral, V., Burton, P., Danesh, J., Downey, P., Elliott, P., Green, J., Landray, M., Liu, B., Matthews, P., Ong, G., Pell, J., Silman, A., Young, A., Sprosen, T., Peakman, T., Collins, R., 2015. UK biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age. PLoS Med 12, e1001779. https://doi.org/10.1371/journal.pmed.1001779 36

Teillac, A., Beaujoin, J., Poupon, F., Mangin, J.-F., Poupon, C., 2017. A Novel AnatomicallyConstrained Global Tractography Approach to Monitor Sharp Turns in Gyri, in: Descoteaux, M., Maier-Hein, L., Franz, A., Jannin, P., Collins, D.L., Duchesne, S. (Eds.), Medical Image Computing and Computer Assisted Intervention – MICCAI 2017. Springer International Publishing, Cham, pp. 532–539. Van Essen, D.C., Smith, S.M., Barch, D.M., Behrens, T.E.J., Yacoub, E., Ugurbil, K., 2013. The WU-Minn Human Connectome Project: An overview. NeuroImage, Mapping the Connectome 80, 62–79. https://doi.org/10.1016/j.neuroimage.2013.05.041 Wang, H., Magnain, C., Wang, R., Dubb, J., Varjabedian, A., Tirrell, L.S., Stevens, A., Augustinack, J.C., Konukoglu, E., Aganj, I., Frosch, M.P., Schmahmann, J.D., Fischl, B., Boas, D.A., 2018. as-PSOCT: Volumetric microscopic imaging of human brain architecture and connectivity. NeuroImage 165, 56–68. https://doi.org/10.1016/j.neuroimage.2017.10.012 Wu, H., Chen, G., Jin, Y., Shen, D., Yap, P.-T., 2016. Embarrassingly Parallel Acceleration of Global Tractography via Dynamic Domain Partitioning. Frontiers in Neuroinformatics 10.

37

Related documents

• Sex-related variability of white-matter tracts is robust to tractography methodology#1038716
• Sex-related variability of white-matter tracts is robust to tractography methodology.#1039955
• fMRI-Targeted High-Angular Resolution Diffusion MR Tractography to Identify Functional Language Tracts in Healthy Controls and Glioma Patients#491782
• What's a White Matter ‘Fiber’?#1090668
• Developmental Trajectory of Brain Structural Connectivity and Cognitive Function From Childhood to Adulthood#847041