Caio Seguin
email: caio.seguin@unimelb.edu.au

Maciej Jedynak

Olivier David

Sina Mansour

Olaf Sporns

Andrew Zalesky

Communication dynamics in the human connectome shape the cortex-wide propagation of direct electrical stimulation

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

Communication between neural elements is paramount to brain function. The complex network of gray matter regions interlinked by white matter fibers-the structural connectome-provides the anatomical substrate for macroscale communication in the brain [START_REF] Sporns | The human connectome: A structural description of the human brain[END_REF][START_REF] Hagmann | Mapping the structural core of human cerebral cortex[END_REF]. Communication via white matter connections underpins interactions between gray matter regions, facilitating the synchronization of neural activity (i.e., functional connectivity) along multiple spatiotemporal scales and giving rise to the brain's rich functional dynamics [START_REF] Honey | Predicting human resting-state functional connectivity from structural connectivity[END_REF][START_REF] Deco | Emerging concepts for the dynamical organization of resting-state activity in the brain[END_REF][START_REF] Park | Structural and functional brain networks: from connections to cognition[END_REF]. Communication processes therefore constitute a bridge between structural and functional descriptions of nervous systems [START_REF] Laughlin | Communication in neuronal networks[END_REF][START_REF] Suárez | Linking structure and function in macroscale brain networks[END_REF].

It is well established that the functional connectivity between two anatomically connected regions is correlated to their structural connectivity [START_REF] Damoiseaux | Greater than the sum of its parts: a review of studies combining structural connectivity and resting-state functional connectivity[END_REF][START_REF] Hermundstad | Structural foundations of resting-state and task-based functional connectivity in the human brain[END_REF][START_REF] Sorrentino | The structural connectome constrains fast brain dynamics[END_REF]. Less well understood, however, are observations of co-activation between anatomically unconnected regions. For example, functionally coupled areas within the same intrinsic restingstate network do not always share a direct structural connection [START_REF] Honey | Can structure predict function in the human brain?[END_REF], as is the case for parts of the parietal cortex and the precuneus involved in the default mode network [START_REF] Van Den Heuvel | Functionally linked resting-state networks reflect the underlying structural connectivity architecture of the human brain[END_REF][START_REF] Alves | An improved neuroanatomical model of the default-mode network reconciles previous neuroimaging and neuropathological findings[END_REF]. Functional connectivity has also been observed between homotopic gray matter loci-normally connected by inter-hemispheric callosal fibers-following complete severance of the corpus callosum [START_REF] Uddin | Residual functional connectivity in the split-brain revealed with resting-state functional mri[END_REF][START_REF] O'reilly | Causal effect of disconnection lesions on interhemispheric functional connectivity in rhesus monkeys[END_REF]. Similarly, targeted pharmacogenetic deactivation of single gray matter structures has been reported to cause global changes in inter-areal co-activation patterns, impacting the activity of regions that were not anatomically connected to the manipulated site [START_REF] Grayson | The rhesus monkey connectome predicts disrupted functional networks resulting from pharmacogenetic inactivation of the amygdala[END_REF]. These examples illustrate that large-scale brain function is not supported exclusively by interactions between directly connected regions, but also by communication dynamics that facilitate polysynaptic information integration between distant and anatomically unconnected regions [START_REF] Suárez | Linking structure and function in macroscale brain networks[END_REF].

Numerous network communication models have been proposed to describe polysynaptic signal propagation in the connectome [START_REF] Avena-Koenigsberger | Communication dynamics in complex brain networks[END_REF]. Computed on structural connectivity, these models quantify the ease of communication between pairs of regions by considering a putative strategy to guide signal transmission. Network communication models provide tractable interpretations of the interplay between structural connectivity and inter-areal interactions, and can thus be used to form and test hypotheses about the relationship between brain structure and function. An emerging body of evidence indicates that these models can explain inter-individual variation in cognitive [START_REF] Seguin | Network communication models improve the behavioral and functional predictive utility of the human structural connectome[END_REF][START_REF] Imms | Navigating the link between processing speed and network communication in the human brain[END_REF] and clinical [START_REF] Raj | A network diffusion model of disease progression in dementia[END_REF][START_REF] Kuceyeski | Longitudinal increases in structural connectome segregation and functional connectome integration are associated with better recovery after mild tbi[END_REF][START_REF] Wang | Synchronization lag in post stroke: relation to motor function and structural connectivity[END_REF][START_REF] De Lange | Shared vulnerability for connectome alterations across psychiatric and neurological brain disorders[END_REF][START_REF] Lella | Communicability distance reveals hidden patterns of alzheimer's disease[END_REF] variables, as well as various aspects of functional and effective connectivity derived from blood-level-oxygen dependent (BOLD) fMRI time courses [START_REF] Grayson | The rhesus monkey connectome predicts disrupted functional networks resulting from pharmacogenetic inactivation of the amygdala[END_REF][START_REF] Goñi | Resting-brain functional connectivity predicted by analytic measures of network communication[END_REF][START_REF] Abdelnour | Network diffusion accurately models the relationship between structural and functional brain connectivity networks[END_REF][START_REF] Abdelnour | Functional brain connectivity is predictable from anatomic network's laplacian eigen-structure[END_REF][START_REF] Osmanlıoglu | System-level matching of structural and functional connectomes in the human brain[END_REF][START_REF] Seguin | Inferring neural signalling directionality from undirected structural connectomes[END_REF][START_REF] Seguin | Network communication models narrow the gap between the modular organization of structural and functional brain networks[END_REF][START_REF] Vázquez-Rodríguez | Gradients of structure-function tethering across neocortex[END_REF][START_REF] Tewarie | Mapping functional brain networks from the structural connectome: Relating the series expansion and eigenmode approaches[END_REF][START_REF] Liu | Time-resolved structure-function coupling in brain networks[END_REF][START_REF] Zamani Esfahlani | Local structure-function relationships in human brain networks across the lifespan[END_REF].

Despite this progress, direct evidence that network communication models reflect aspects of biological neural signaling remains lacking. Slow-fluctuating BOLD time series-the focus of the vast majority of studies to date-are unable to capture inter-areal communication the at millisecond time scales inherent to signal propagation. Conversely, high-temporal-resolution electrophysiological recordings are typically limited in their ability to measure spatially-resolved neural activity across the whole brain [START_REF] Palva | Ghost interactions in meg/eeg source space: A note of caution on inter-areal coupling measures[END_REF][START_REF] Parvizi | Promises and limitations of human intracranial electroencephalography[END_REF]. These technological obstacles have thus far hampered the biological validation of brain network communication models.

Here, we address this challenge using invasive stereotactic electroencephalography (SEEG) recordings of cortico-cortical evoked potentials (CCEPs) following direct electrical stimulation of the human brain. We use data from the F-TRACT project [START_REF] David | Probabilistic functional tractography of the human cortex[END_REF][START_REF] Trebaul | Probabilistic functional tractography of the human cortex revisited[END_REF][START_REF] Lemaréchal | A brain atlas of axonal and synaptic delays based on modelling of corticocortical evoked potentials[END_REF] to derive wholebrain, high-spatial-resolution maps of stimulus propagation measured at millisecond resolution. We hypothesize that modeling network communication on the human structural connectome can predict the cortex-wide propagation of focal electrical stimulation. Importantly, we consider a wide range of models-from efficient routing protocols to passive diffusion processes-to determine which conceptualizations of network communication best explain stimulus transmission. In parallel, we leverage the interpretability of our modeling framework to make inferences about which organizational properties of the connectome contribute to inter-areal communication. We posit that progress in understanding the spread of exogenous, interventional stimulation will provide insights into the underlying mechanisms governing large-scale neural communication.

RESULTS

The F-TRACT project [START_REF] David | Probabilistic functional tractography of the human cortex[END_REF][START_REF] Trebaul | Probabilistic functional tractography of the human cortex revisited[END_REF] processed multicenter SEEG data acquired in a large sample of patients with drug-resistant epilepsy. The position and number of implanted electrodes were separately determined for each patient based on prior clinical knowledge. Patients underwent focal single-pulse stimulation and invasive recordings of neural activity at millisecond resolution. Compiling data from a total of 550 participants, 29,055 stimulations and 2.77 million pairs of recording electrodes, allowed the computation of whole-brain, group-level probability (P) and amplitude (A) matrices of regional activation following stimulation (Fig 1a; Materials and Methods). The matrix entries P ij ∈ [0, 1] and A ij ∈ R >0 capture, respectively, the probability of observing a significant excitatory response in region j following stimulation of region i, and the median amplitude of significant responses.

To investigate the role of brain connectivity on stimulus response, we mapped structural and functional brain networks from MRI data of 1,000 healthy participants of the Human Connectome Project [START_REF] Van Essen | The WU-Minn human connectome project: an overview[END_REF] (HCP; Fig 1b; Materials and Methods). We generated four grouprepresentative normative measures: structural connectivity (SC), structural connection frequency (F) across HCP participants, functional connectivity (FC) and interregional Euclidean distance (ED). While our primary goal is to understand communication via white matter connectivity, considering functional and spatial factors enabled a comprehensive examination of how stimulus propagation is associated with multiple aspects of brain organization.

Whole-brain maps of electrical stimulus propagation

We analyzed response probability estimates between 4,989 region pairs (15.4% of all intra-hemispheric region pairs; Fig 1c), out of which 2,792 were likely to be anatomically connected (F > 75%; violet matrix entries) and 1,065 were likely to be anatomically unconnected (F < 25%; copper matrix entries). Stimulation and recording sites covered the majority of the cortical surface (Fig 1d; 81% and 84% of regional coverage, respectively), and were primarily concentrated in areas involved in refractory focal epilepsy. For amplitude, we analyzed estimates between 857 region pairs (2.6% of all intra-hemispheric region pairs; Fig 1g), out of which 644 and 96 were likely and unlikely to be anatomically connected. Stimulation and recording sites remained distributed across large portions of the cortex (Fig 1h; 47% and 54% of regional coverage, respectively). Note that amplitude measurements are sparser since they only encode stimulations that elicited significant responses. Despite the heterogeneous spatial distribution of electrode implantations, pairwise response measures were distributed across multiple regions and systems, thus allowing for whole-brain analyses of stimulus propagation at high spatiotemporal resolution.

We found that anatomically connected regions showed higher response probability (Wilcoxon rank sum test p < 10 ,j, probability and amplitude were correlated to the SC (e.g., ρ con = 0.46, 0.36 for P and A, respectively), FC (e.g., ρ all = 0.42, 0.34), and especially the ED (e.g., ρ all = -0.62, -0.50) between stimulation and recording sites (all p < 10 -32 ).

In line with previous literature [START_REF] Betzel | Structural, geometric and genetic factors predict interregional brain connectivity patterns probed by electrocorticography[END_REF], these results indicate that multiple aspects of connectome architecture may shape stimulus propagation, including anatomical connectivity [START_REF] Sorrentino | The structural connectome constrains fast brain dynamics[END_REF][START_REF] Stiso | White matter network architecture guides direct electrical stimulation through optimal state transitions[END_REF], resting-state functional coupling [START_REF] Veit | Temporal order of signal propagation within and across intrinsic brain networks[END_REF][START_REF] Keller | Intrinsic functional architecture predicts electrically evoked responses in the human brain[END_REF], and brain geometry [START_REF] Trebaul | Probabilistic functional tractography of the human cortex revisited[END_REF].

Modeling brain network communication

We next sought to investigate the network communication dynamics that facilitate polysynaptic (i.e., multistep) stimulus propagation (Fig 2a). A network communication model determines a policy to guide signals through structural connectivity, and can therefore be used to quantify properties of stimulus transmission between regions. We considered five popular measures of brain network communication: shortest path efficiency (a) Schematics of the F-TRACT project. Subjects were implanted with SEEG electrodes capable of delivering localized stimuli and recording electrical activity at millisecond resolution. The toy example shows a stimulation performed on a single subject, where the effects of a single-pulse stimulus are recorded (schematics are illustrative and do not reflect the bipolar stimulation montage and the true spacing between contacts). A significant response is measured by contact 1, as defined by a z-scored CCEP larger than 5. Contact 2 does not record a significant response. Based on its spatial coordinate, each contact was assigned to a cortical gray matter region. Whole-brain matrices of response probability (P) and median amplitude (A) were computed by combining data from 29,055 stimulations in 550 subjects. (b) Normative, group-representative measures of brain connectivity were inferred from 1,000 healthy participants of the Human Connectome Project using structural, diffusion, and resting-state functional MRI data. (c) Pairs of regions used to model response probability. Violet and copper matrix entries denote regions that are likely and unlikely to be anatomically connected, respectively. Warm-and cool-colored line plots indicate the number of times each region featured as a stimulation and recording site, respectively, and were projected onto the cortical surface on panel (d). (e) Response probability distributions for all (black), connected (violet), and unconnected (copper) region pairs. (f ) Bars show the absolute value of Spearman's ρ between response probability and ED, SC and FC, computed for all (left), connected (center), and unconnected (right) region pairs. Circles and error bars indicate the mean and 95% confidence interval correlations obtained from 1,000 bootstrapped samples. (g-j) Same as (c-f) but for response amplitude. [START_REF] Latora | Efficient behavior of smallworld networks[END_REF][START_REF] Kaiser | Nonoptimal component placement, but short processing paths, due to longdistance projections in neural systems[END_REF], navigation efficiency [START_REF] Boguna | Navigability of complex networks[END_REF][START_REF] Seguin | Navigation of brain networks[END_REF], search information [START_REF] Goñi | Resting-brain functional connectivity predicted by analytic measures of network communication[END_REF][START_REF] Rosvall | Searchability of networks[END_REF], communicability [START_REF] Andreotti | Validation of network communicability metrics for the analysis of brain structural networks[END_REF][START_REF] Estrada | Communicability in complex networks[END_REF], and diffusion efficiency [START_REF] Goñi | Exploring the morphospace of communication efficiency in complex networks[END_REF] (see Fig 2b for brief descriptions of measures and the Materials and Methods for technical details). Covering a wide range of signaling conceptualizations, these measures are positioned along a spectrum from routingcommunication via efficient, selectively accessed pathsto diffusion-spreading along multiple network fronts [START_REF] Avena-Koenigsberger | A spectrum of routing strategies for brain networks[END_REF].

Importantly, models located at different positions of the routing-diffusion spectrum rely on distinct properties of the underlying structural connectivity to guide signaling (Fig 2c). Routing measures (shortest path efficiency and navigation efficiency) consider that stimuli are communicated exclusively via efficient paths comprising a small number of high-weight connections. In contrast, diffusion measures (communicability and diffusion efficiency) posit that the broader network topology also plays a significant role in communication. For example, the presence of multiple alternative paths from stimulation to recording sites is expected to increase the strength of stimulus response. Lastly, search information combines aspects of both routing and diffusion. In the context of brain stimulation, it can be interpreted as the propensity for a stimulus to be dispersed and diluted while traveling via the shortest path from stimulation to recording sites.

Determining whether certain measures yield better explanations of stimulus response than others can therefore provide insight into which topological properties of the connectome are most relevant for neural communication.

Network communication models explain polysynaptic propagation of electrical stimulation

We correlated network communication measures computed on the structural connectome with response probability and amplitude. To account for the contribution of spatial proximity between stimulation and recording sites, we computed both full (ρ) and partial (ρ ) correlations, the latter controlling for the effect of ED.

We first observed that, despite the strong influence of ED, connectivity-based measures explained significant variance in stimulus propagation that was not accounted for by spatial proximity ( cantly greater variability in stimulus propagation than structural connectivity alone (Fig S1). In particular, search information consistently ranked as the most explanatory network measure of response probability and amplitude (e.g., ρ all = -0.60, -0.47, Figs 2c,h, respectively; both p < 10 -32 ), in many cases yielding statistically greater correlations than alternative measures (Fig S1). Interestingly, top-performing network measures, such as search information and communicability, surpassed the explanatory power of SC weights even for associations computed exclusively for pairs of anatomically connected regions (Fig S1b,e). This suggests that stimulus transmission between connected sites may not be supported entirely by direct fibers, but also facilitated by propagation through intermediate regions and connections in the surrounding network topology.

Response amplitude measurements between anatomically unconnected regions describe instances of significant multi-step stimulus transmission through the connectome, thus reflecting a primary conceptual application of network communication models (Fig 2a,b). Crucially, in this scenario (Fig 3f), search information and communicability markedly outperformed alternative network measures and ED (SI: ρ unc = -0.48, p = 8×10 -7 , CMY: ρ unc = 0.42, p = 2 × 10 -5 ; Fig S1f,l). These results provide evidence that network communication measures can explain the propagation of electrical stimuli-traveling at millisecond resolution-between anatomically unconnected regions.

Multivariate prediction of electrical stimulus propagation

Next, we combined SC, FC, ED, and all network communication measures into a single multivariate predictive model (Materials and Methods). Using a machine learning method that penalizes model complexity, we sought to identify sparse sets of predictors that explained nonoverlapping variance in response probability and amplitude. We found that a simple predictive model comprising Euclidean distance (40.2% average contribution), search information (32.5%) and FC (24.6%) explained R 2 = 54% of out-of-sample variance in response probability (Fig 4a-c). For amplitude, a combination of Euclidean distance (56.5% average contribution), communicability (29.5%), as well as search information and FC (14% joint contribution), yielded out-of-sample R 2 = 31% (Fig 4d-f).

These results corroborate and complement our univariate analyses, indicating that a parsimonious combination of structural, functional and spatial factors can accurately predict both the likelihood and the strength of cortex-wide responses to focal brain stimulation. Additionally, these analyses further indicate that search information and communicability provide the most explanatory conceptualizations of stimulus propagation through structural connectivity.

Control analyses

We performed multiple control analyses (Materials and Methods) and found that our results remained consistent across a number of methodological choices (Fig S2), including replication in an independently acquired and processed normative connectivity dataset (Figs S3, S4). In addition, using ensembles of randomized null networks, we observed that our findings were contingent on the empirical topology of SC, such that rewiring of connections markedly reduced correlations between network communication measures and stimulus response (Fig S5).

DISCUSSION

Here, we used a large dataset of direct electrical stimulation to study inter-areal communication in the human brain. Considering the causal effects of interventional stimulation allowed us to characterize the propagation of electrical pulses without the need for statistical or information-theoretical measures of information transfer between regions. In addition, we leveraged intraand inter-subject variability in electrode placement to derive whole-brain, high-spatiotemporal-resolution maps of stimulus propagation. Together, these factors enabled us to test whether computational modeling of communication dynamics in the human connectome can explain empirical stimulus transmission through the brain.

Towards the biological validation of brain network communication models

Taken together, our analyses identified two key findings. First, we provide evidence that network communication models-computed on structural connectivity in-ferred from tractography and MRI-reflect patterns of spatially-resolved, millisecond resolution, stimulus propagation measured across the whole brain. Our results recapitulate previous findings that structural connectivity shapes fast brain dynamics [START_REF] Sorrentino | The structural connectome constrains fast brain dynamics[END_REF], and corroborate the use of network communication models to study neural information processing in healthy [START_REF] Mišić | Cooperative and competitive spreading dynamics on the human connectome[END_REF][START_REF] Vézquez-Rodríguez | Signal propagation via cortical hierarchies[END_REF] and clinical populations [START_REF] Raj | A network diffusion model of disease progression in dementia[END_REF][START_REF] Kuceyeski | Longitudinal increases in structural connectome segregation and functional connectome integration are associated with better recovery after mild tbi[END_REF][START_REF] Wang | Synchronization lag in post stroke: relation to motor function and structural connectivity[END_REF][START_REF] De Lange | Shared vulnerability for connectome alterations across psychiatric and neurological brain disorders[END_REF][START_REF] Lella | Communicability distance reveals hidden patterns of alzheimer's disease[END_REF].

Second, adding to previous literature [START_REF] Grayson | The rhesus monkey connectome predicts disrupted functional networks resulting from pharmacogenetic inactivation of the amygdala[END_REF][START_REF] Goñi | Resting-brain functional connectivity predicted by analytic measures of network communication[END_REF][START_REF] Betzel | Structural, geometric and genetic factors predict interregional brain connectivity patterns probed by electrocorticography[END_REF][START_REF] Andreotti | Validation of network communicability metrics for the analysis of brain structural networks[END_REF], we found that-across multiple control analyses-search information and communicability led to the most accurate and robust accounts of stimulus propagation. This included correlations to the response amplitude between anatomically unconnected regions, indicating that these measures can explain multi-step propagation through the connectome. Despite being computed on the same structural substrate, alternative measures yielded weaker or inconsistent associations, and did not provide additional explanatory power when considered in multivariate statistical models. This specificity in predictive utility constitutes an important step towards determining which network communication models are more suitable to describe biological neural signaling.

Our work complements parallel research on network control theory [START_REF] Muldoon | Stimulationbased control of dynamic brain networks[END_REF][START_REF] Stiso | White matter network architecture guides direct electrical stimulation through optimal state transitions[END_REF] and biophysical models of neural activity [START_REF] Gollo | Mapping how local perturbations influence systems-level brain dynamics[END_REF][START_REF] Cabral | Functional connectivity dynamically evolves on multiple time-scales over a static structural connectome: Models and mechanisms[END_REF][START_REF] Papadopoulos | Relations between large-scale brain connectivity and effects of regional stimulation depend on collective dynamical state[END_REF]. Recent reports have highlighted similarities and opportunities for synergy between these approaches for modeling brain function [START_REF] Palmigiano | Flexible information routing by transient synchrony[END_REF][START_REF] Srivastava | Models of communication and control for brain networks: distinctions, convergence, and future outlook[END_REF]. As such, the present evidence in support of network communication models also bolsters research on these alternative directions. In particular, our work calls for the extension of efforts towards the empirical and biological validation of computational models of macroscale brain function.

Understanding how brain network organization facilitates inter-areal stimulus propagation

A key feature of network communication models is that they provide tractable interpretations of the interplay between brain structure and function. Our findings reveal insights into how topological properties of structural connectivity facilitate inter-areal signaling (Fig 2). Search information and communicability conceptualize network communication as a diffusive process, in which stimulus propagation is shaped by signal dispersion and transmission via multiple network fronts. This notion challenges the mainstay assumption in network neuroscience that communication occurs exclusively via shortest paths [START_REF] Seguin | Navigation of brain networks[END_REF][START_REF] Goñi | Exploring the morphospace of communication efficiency in complex networks[END_REF], and suggests that polysynaptic dynamics may play a role even for anatomically connected regions.

We envision that these results may contribute to the design of brain stimulation paradigms that-beyond maps of direct connectivity-also leverage knowledge of network communication dynamics [START_REF] Horn | Opportunities of connectomic neuromodulation[END_REF]. For example, models of polysynaptic transmission could inform the selection of stimulation sites for indirect modulation of deep gray matter structures by non-invasive techniques, such as transcranial magnetic stimulation [START_REF] Cash | Using brain imaging to improve spatial targeting of transcranial magnetic stimulation for depression[END_REF][START_REF] Sydnor | Cortical-subcortical structural connections support transcranial magnetic stimulation engagement of the amygdala[END_REF].

In addition to propagation constrained by white matter, spatial and functional factors were also predictive of stimulus response. We note that while the influence of spatial proximity may reflect neurophysiological evidence that amplitude decays as a function of axonal length [START_REF] Trebaul | Probabilistic functional tractography of the human cortex revisited[END_REF], it could also be confounded by factors such as non-axonal conduction, spatial autocorrelations in MRI data [START_REF] Poldrack | Handbook of functional MRI data analysis[END_REF], and inaccurate reconstructions of short-range fibers by tractography [START_REF] Zhang | Quantitative mapping of the brain's structural connectivity using diffusion MRI tractography: A review[END_REF]. We also found that, despite a relatively weaker explanatory power, functional connectivity explained unique variance in stimulus propagation. In line with previous studies [START_REF] Veit | Temporal order of signal propagation within and across intrinsic brain networks[END_REF][START_REF] Keller | Intrinsic functional architecture predicts electrically evoked responses in the human brain[END_REF], this indicates that intrinsic functional dynamics may relate to stimulus response in ways complementary to structural and spatial factors [START_REF] Papadopoulos | Relations between large-scale brain connectivity and effects of regional stimulation depend on collective dynamical state[END_REF]. Multimodal models of connectome communication that incorporate information on both structural and functional connectivity may therefore provide a fruitful direction for future research [START_REF] Zheng | Local vulnerability and global connectivity jointly shape neurodegenerative disease propagation[END_REF][START_REF] Amico | Toward an information theoretical description of communication in brain networks[END_REF].

Limitations and future directions

A number of limitations of the presence work should be discussed. First, stimulus response were measured in individuals with refractory epilepsy. While conservative methods were applied to filter out pathological activity from our analyses (see Materials and Methods), the extent to which recordings from epilepsy patients are representative of healthy brain activity remains an active topic of research [START_REF] Parvizi | Promises and limitations of human intracranial electroencephalography[END_REF][START_REF] Lemaréchal | A brain atlas of axonal and synaptic delays based on modelling of corticocortical evoked potentials[END_REF]. Similarly, CCEPs elicited by direct electrical stimulation may not reflect physiological patterns of inter-areal communication. While future work is required to investigate these issues, we note that our findings are congruent with previous reports that diffusion-based communication models explained activity measured using electrocorticography in the absence of electrical stimulation [START_REF] Betzel | Structural, geometric and genetic factors predict interregional brain connectivity patterns probed by electrocorticography[END_REF] and non-invasive magnetoencephalography in healthy adults [START_REF] Tewarie | Predicting timeresolved electrophysiological brain networks from structural eigenmodes[END_REF].

Our analyses focused on intra-hemispheric stimulus propagation (Fig. 1c,g). This was a result of selecting region pairs for which stimulus response measures were reliably estimated based on a large number of stimulations (see Materials and Methods). Nonetheless, determining whether our results extend to inter-hemispheric communication is an important step for subsequent research. Lastly, while our results corroborate the use of grouplevel normative connectivity to study neuromodulation in clinical populations [START_REF] Horn | Connectivity predicts deep brain stimulation outcome in parkinson disease[END_REF], future research on modeling subject-level data is necessary.

Conclusion

In conclusion, through the analyses of more than 2.7 million SEEG recordings, we provide novel insights into connectome communication and polysynaptic stimulus propagation. This work is a step towards the development and validation of biologically realistic models of brain communication. Our results provide new avenues for research on macroscale neural information processing and brain stimulation paradigms that leverage network communication dynamics.

MATERIALS AND METHODS

F-TRACT dataset

The F-TRACT project is a continuing effort in collecting, processing and compiling iEEG recordings of drug-resistant epilepsy patients being prepared for resective surgery [START_REF] David | Probabilistic functional tractography of the human cortex[END_REF][START_REF] Trebaul | Probabilistic functional tractography of the human cortex revisited[END_REF]. In this study, we used data from 550 patients (50% women, age 31 ± 10 years old) collected across 20 epilepsy centers in Europe, North America and Asia (https://f-tract.eu/consortium). Patients were implanted with intracerebral depth electrodes (stereo-electroencephalography; SEEG) with an average of 87 ± 37 electrode contacts per participant. Corticocortical evoked potentials (CCEPs) were recorded following low-frequency, bipolar direct electrical stimulation with the following parameters: 1 Hz stimulation, biphasic (79% of CCEPs) or monophasic (21%) electrical pulses, 3.3 ± 1.2 mA mean pulse intensity, 1.0 ± 0.4 ms mean pulse duration, and 3.3 ± 1.5 µC mean pulse charge.

CCEPs were processed as follows. First, SEEG signals were preprocessed to identify stimulation events, filter bad channels [START_REF] Tuyisenge | Automatic bad channel detection in intracranial electroencephalographic recordings using ensemble machine learning[END_REF], and correct for stimulation artifacts [73]. Recordings with high likelihood of reflecting pathological activity were identified and discarded using Delphos software to mitigate the impact of epileptogenic processes in further analyses [START_REF] Lemaréchal | A brain atlas of axonal and synaptic delays based on modelling of corticocortical evoked potentials[END_REF][START_REF] Roehri | Time-Frequency strategies for increasing oscillation detectability in intracerebral[END_REF]. CCEPs on each electrode contact were baseline corrected and z-scored based on the [-400, -10] ms time window before stimulation. A stimulus response was deemed significant if a z-CCEP above z = 5 was observed within 800 ms of the stimulation. The response amplitude of a single stimulation was defined as the first z-CCEP peak above the significance threshold.

The position of electrode contacts was computed based on a patient's structural MRI (T1-weighted imaging) and computerized tomography scans. Contacts placed in gray matter were assigned to the corresponding regions of the HCP-MMP1.0 parcellation [START_REF] Glasser | A multi-modal parcellation of human cerebral cortex[END_REF]. For contacts placed in white matter, we first computed a 3 mm radius sphere centered in the contact's coordinates and assigned them to the gray matter parcel with highest voxel representation within the sphere.

Whole-brain, group-level response probability (P) and amplitude (A) matrices were computed by compiling CCEP data from 2,769,985 SEEG recordings following 29,055 electrical stimulations across 550 individuals. The matrix entries P ij ∈ [0, 1] and A ij ∈ R >0 capture, respectively, the probability of observing a significant excitatory response in region j following stimulation of region i, and the median z-scored CCEP amplitude of significant Note that P and A are asymmetric, such that P ij may differ to P ji . Data from the F-TRACT project is openly available at https://f-tract.eu/ and further details on data acquisition and processing are provided in [START_REF] David | Probabilistic functional tractography of the human cortex[END_REF][START_REF] Trebaul | Probabilistic functional tractography of the human cortex revisited[END_REF][START_REF] Lemaréchal | A brain atlas of axonal and synaptic delays based on modelling of corticocortical evoked potentials[END_REF].

Normative measures of brain connectivity: HCP dataset

Structural connectivity (SC) was inferred from minimally preprocessed, high-resolution diffusion-weighted MRI data from 1,000 healthy young adults participating in the Human Connectome Project. Details of diffusion MRI acquisition and preprocessing are described in [START_REF] Sotiropoulos | Advances in diffusion MRI acquisition and processing in the human connectome project[END_REF][START_REF] Glasser | The minimal preprocessing pipelines for the human connectome project[END_REF]. Whole-brain white matter tractograms of individual participants were mapped using a probabilistic tractography pipeline (MRtrix3 software [START_REF] Tournier | MRtrix: Diffusion tractography in crossing fiber regions[END_REF], multi-shell multi-tissue constrained spherical deconvolution [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[END_REF], anatomically constrained tractography [START_REF] Smith | Anatomically-constrained tractography: im-proved diffusion MRI streamlines tractography through effective use of anatomical information[END_REF], iFOD2 tracking algorithm [START_REF] Tournier | Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions[END_REF], 5M streamlines seeded from the gray-white matter interface; see [START_REF] Mansour | High-resolution connectomic fingerprints: Mapping neural identity and behavior[END_REF] for further details). Cortical gray matter was divided into 360 regions using the HCP-MMP1.0 parcellation [START_REF] Glasser | A multi-modal parcellation of human cerebral cortex[END_REF], while the Freesurfer parcellation [START_REF] Fischl | Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain[END_REF] was used to define 14 subcortical structures (left and right thalamus, caudate, putamen, pallidum, hippocampus, amygdala and accumbens area). Note that while stimulus response measures were available only in the cortex, the inclusion of subcortical regions in structural connectomes enables the modeling of propagation mediated by the subcortex. The structural connectivity weight between two gray matter regions was defined as the number of streamlines between them divided by the sum of their surface areas, approximated as the number of voxels in the gray-white matter inter-face of each region [START_REF] Seguin | Inferring neural signalling directionality from undirected structural connectomes[END_REF]. Connections comprising fewer than 5 streamlines were pruned to reduce the high false positive rate of probabilistic tractography [START_REF] Zalesky | Connectome sensitivity or specificity: which is more important?[END_REF]. A grouprepresentative connectivity matrix was computed using a consensus method that combined individual-level networks while preserving the average connection density across subjects [START_REF] Betzel | Distance-dependent consensus thresholds for generating group-representative structural brain networks[END_REF]. This resulted in a final 374 × 374 SC matrix with 25% connection density.

Functional connectivity (FC) was mapped using minimally preprocessed, ICA-FIX resting-state functional MRI data from the same sample of 1,000 HCP participants. For each subject, four 14 minutes and 33 second scans (0.72 second TR) were acquired following details described in [START_REF] Glasser | The minimal preprocessing pipelines for the human connectome project[END_REF][START_REF] Smith | Resting-state fMRI in the human connectome project[END_REF]. Voxel-level blood-oxygen-leveldependent time series were linearly detrended, band-pass filtered and standardized [START_REF] Parkes | An evaluation of the efficacy, reliability, and sensitivity of motion correction strategies for resting-state functional MRI[END_REF]. Following previous work [START_REF] Satterthwaite | An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data[END_REF], a total of four nuisance variables were regressed out: the global signal (GS), the GS squared, the GS derivative, and the squared GS derivative. The time series of voxels located in the same gray matter parcel were averaged and FC was computed as the Pearson correlation between regional time series. A final 374 × 374 grouprepresentative FC matrix was computed as the average of 4,000 individual matrices (4 sessions per subject for 1,000 subjects).

Normative measures of brain connectivity: MICA dataset

The Multimodal Imaging and Connectome Analysis (MICA) dataset is an openly available resource (https: //portal.conp.ca/) containing structural and functional brain networks for 50 healthy adults (23 women; 29.5 ± 5.6 years old) [START_REF] Royer | An open MRI dataset for multiscale neuroscience[END_REF]. Brain networks were inferred using independent acquisition, preprocessing and connectivity mapping pipelines. We utilized the MICA brain networks comprising 374 regions, delineated according to the same cortical and subcortical parcellations used in the HCP dataset. To account for the nearly fullyconnected nature of MICA structural networks (average connection density 84.1% across subjects) and prune potential false positives [START_REF] Zalesky | Connectome sensitivity or specificity: which is more important?[END_REF], we thresholded individual-level SC matrices by keeping only the top-25% strongest structural connections. Group-representative structural and functional networks were then derived following the same procedures described above for the HCP dataset.

Network communication models

Let W ∈ R N ×N denote the matrix of SC weights between N regions, where W ij provides a measure of the strength and reliability of the anatomical connection between regions i and j. Anatomically connected and unconnected region pairs have W ij > 0 and W ij = 0, respectively. We define the matrix of SC lengths L = 1/W , where L ij quantifies the distance of travel cost between i and j. The remapping of connection weights into lengths is required for the computation of network communication models that aim to minimize the cost of transmitting signals between regions. Measures of network communication were computed using the Brain Connectivity Toolbox [START_REF] Rubinov | Weight-conserving characterization of complex functional brain networks[END_REF] and were defined as follows.

Shortest path efficiency

The Floyd-Warshal algorithm was used to identify the sequence of regions {i, u, ..., v, j} such that Λ * ij = L iu +. . . +L vj is the minimum transmission cost between regions i and j. Shortest path efficiency was defined as [START_REF] Latora | Efficient behavior of smallworld networks[END_REF].

SP E ij = 1/Λ * ij

Navigation efficiency

Navigation is a greedy routing protocol that propagates signals based on a measure of inter-regional distance [START_REF] Boguna | Navigability of complex networks[END_REF]. Following previous work [START_REF] Seguin | Navigation of brain networks[END_REF], navigation from a source region i to a target region j was implemented as follows. First, identify i's neighbor with the shortest Euclidean distance to j and progress to it. Repeat this process until j is reached (successful navigation) or a region is revisited (failed navigation). The length of successful navigation paths is given by Λ ij = L iu +. . . +L vj , where {i, u, ..., v, j} is the sequence of nodes visited during routing, whereas failed navigation entails Λ ij = ∞. Navigation efficiency was defined as N E ij = 1/Λ ij .

Diffusion efficiency

The diffusion efficiency quantifies the ease of communication from region i to region j under the assumption that signals are transmitted via unbiased random walks [START_REF] Goñi | Exploring the morphospace of communication efficiency in complex networks[END_REF]. In this scenario, the probability of traveling from p to q in one step is given by the transition matrix T pq = W pq / N n=1 W pi . The transition probability matrix can be used to analytically compute H ij , the average number of steps required for a random walker to travel from i to j [START_REF] Zhou | Network landscape from a brownian particle's perspective[END_REF]. Diffusion efficiency is defined as

DE ij = 1/ H ij .

Search information

Search information quantifies the amount of information necessary to bias a random walker to travel via shortest paths, capturing the accessibility of efficient routes under a diffusive model of communication [START_REF] Goñi | Resting-brain functional connectivity predicted by analytic measures of network communication[END_REF][START_REF] Rosvall | Searchability of networks[END_REF]. Considering the transition matrix T of a random walk process, P i ij = T iu ×. . . ×T vj is the probability of a random walker serendipitously traveling from region i to region j via their shortest path {i, u, ..., v, j}. Search information was defined as SI ij = -log 2 (Π ij ).

Communicability

Communicability is defined as a weighted sum of the total number of walks between two regions, in which the contribution of each walk to signal propagation is inversely proportional to their length [START_REF] Estrada | Communicability in complex networks[END_REF][START_REF] Crofts | A weighted communicability measure applied to complex brain networks[END_REF]

. For- mally, CM Y ij = ∞ n=0 W n ij /n! = e W ij , where W ij = W ij /( √ s i √ s j
) is a normalized connectivity matrix that attenuates the influence of high strength region and si = N n=1 s n is the strength of region i [START_REF] Crofts | A weighted communicability measure applied to complex brain networks[END_REF].

Univariate modeling of stimulus response and control analyses SC, FC, ED, and the 5 network communication measures were correlated to response probability and amplitude using Spearman rank correlation. Associations were computed across pairs of cortical regions. In order to reduce potential effects of measurement noise and interindividual variability in CCEPs, we considered, in the main text, pairs of regions for which (i) a minimum of 100 stimulation experiments were available for the computation of probability, (ii) a minimum of 100 significant stimulation responses were available for the computation of amplitude, and (iii) SEEG recordings were across a minimum of 5 different patients. In addition, to mitigate potential confounding factors relating to non-axonal volume conduction, spatial autocorrelations in MRI data [START_REF] Poldrack | Handbook of functional MRI data analysis[END_REF], and inaccurate reconstructions of short-range U-fibers ny tractography [START_REF] Zhang | Quantitative mapping of the brain's structural connectivity using diffusion MRI tractography: A review[END_REF], the results in the main text excluded pairs of regions closer than 20 mm in Euclidean distance [START_REF] Power | Functional network organization of the human brain[END_REF].

We performed a number of control analyses to determine whether our findings were robust to variations in the methodological choices described above. Two sets of control analyses were considered. Firstly, in Fig S2, we systematically analyzed variations to the following parameters:

• Minimum number of stimulation experiments/significant responses between region pairs (100 [main text] vs. 50 [control analysis]).

• Minimum structural connectivity frequency across HCP subjects for considering two regions to be likely anatomically connected (75% [main text] vs. 90% [control analysis]).

• Maximum structural connectivity frequency across HCP subjects for considering two regions unlikely anatomically connected (25% [main text] vs. 10% [control analysis]).

• z-scored CCEP threshold for significant responses (z = 5 [main text] vs. z = 4, 6 [control analyses]).

• Minimum Euclidean distance between regions pairs (20 mm [main text] vs. 0 mm [control analyses]).

Secondly, in Figs S3, we re-examined all the scenarios above starting from structural and functional connectivity from the Multimodal Imaging and Connectome Analysis (MICA) dataset.

Multivariate machine learning predictive model

Let y ∈ R n×1 be a vector of outcome variables corresponding to the response probability or amplitude of n pairs of regions. Predictive modeling was carried out considering all pairs of regions used in the univariate analyses, so that n = 4989 for probability and n = 857 for amplitude. Let X ∈ R n×8 be a matrix of predictor variables comprising the SC, FC, ED, and the 5 measures of network communication between n pairs of regions. We used a lasso regression [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] to predict y from X.

The lasso regression was chosen as it implements a feature selection procedure aimed at identifying a sparse set of variables in X that accurately predicts y. This framework seeks to account for correlations between the variables in X by culling measures that offer redundant predictive utility. As such, predictor variables selected by the lasso regression tend to explain complementary portions of the variance in y.

Model training was done using 10-fold cross-validation by splitting X and y into train and test sets. Model fitting was performed exclusively on train sets and prediction accuracy was evaluated exclusively on test sets. To account for sensitivity to particular train-test splits, we repeated the 10-fold cross validation routine 15 times [START_REF] Li | Global signal regression strengthens association between resting-state functional connectivity and behavior[END_REF]. Computations were carried out using the cvglmnet MATLAB function [START_REF] Qian | Glmnet for matlab[END_REF].

More specifically, let {X a , y a } and {X e , y e } denote a split of X and y into train and test sets, respectively. The lasso regression computes β ∈ R 8×1 as min

β { 1 n ||y a -X a β|| 2 2 + λ||β|| 1 },
where the hyperparameter 0 ≤ λ ≤ 1 determines the sparsity of feature selection. The value of λ was tuned using nested cross-validation, in which each outer training set X a was further split into 10 inner training sets. For each inner training fold, we considered K logarithmically spaced λ from λ min to λ max , where λ min = λ max and λ max is the smallest parameter value such that ||β inner || 1 = 0. Values of K and were dynamically determined by the nested cross-validation algorithm, with typical values around K = 100, and = 0.01 [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF]. The mean squared error (MSE) of the resulting K models were then evaluated on the inner test folds. To promote sparsity, the largest λ with MSE within one standard error of the optimal MSE was selected [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF]. The selected λ was then used to compute β for the training set X a . The prediction accuracy of each outer folder was evaluated as the Pearson correlation and the root mean squared error (RMSE) between ŷe and y e , where ŷe = X e β. The contribution of predictor i was computed as

|β i |/ 8 i |β i |.
This procedure was repeated for 150 pairs of outer train and test sets (15 repetitions of 10-fold cross validation).

Rewiring null models SC topology was randomized using standard and costpreserving versions of the Maslov-Sneppen rewiring algorithm [START_REF] Maslov | Specificity and stability in topology of protein networks[END_REF]. Standard rewiring was performed by swapping pairs of connections while preserving empirical weight distribution and degree sequence. Cost-preserving rewiring was implemented by restricting swaps to pairs of connections with approximately the same length [START_REF] Betzel | Specificity and robustness of long-distance connections in weighted, interareal connectomes[END_REF], as defined by the Euclidean distance between nodes. This routine preserves empirical weight distribution, degree sequence, connection length distribution, and connection weight-length relationship. In both routines, each connection was rewired once on average. Code implementing the standard and cost-preserving routines are available through references [START_REF] Rubinov | Weight-conserving characterization of complex functional brain networks[END_REF] and [START_REF] Betzel | Specificity and robustness of long-distance connections in weighted, interareal connectomes[END_REF], respectively.

An ensemble of 1,000 standard and cost-preserving randomized structural networks was generated from empirical SC. Network communication measures were computed for each randomized network and correlated to response probability and amplitude to form distributions of null correlations (Fig S5). A network communication measure was considered to outperform a null model if its empirical correlation to stimulus response exceeded the 95th percentile of the null correlation distribution (nonparametric significance test at α = 5%). Evidence that a network communication measure outperforms a null model indicates that its association to stimulus response is not trivially explained by the SC attributes preserved during rewiring, and is instead contingent on the empirical topology of anatomical connectivity. Control analysis in which response to stimuli were considered significant for z-scored CCEP amplitude larger than 4. (i-j) Control analysis in which response to stimuli were considered significant for z-scored CCEP amplitude larger than 6. (k-l) Control analysis in which no restrictions were imposed on Euclidean distance between pairs of regions. FIG. S3. Control analyses for normative brain connectivity inferred from the Multimodal and Connectome Analysis (MICA) dataset. Each row of plots shows results obtained for a different set of methodological choices. Panels (a-b) shows results obtained considering the parameters utilized in the main text: (i) pairs of regions were included in the analyses of response probability (amplitude) if they had at least 100 stimulation experiments (significant responses); (ii) pairs of were considered to be anatomically connected (unconnected) if they shared a connection in at least 75% (no more than 25%) of participants; (iii) response to stimuli were considered significant for z-scored CCEP amplitude larger than 5; and (iv) pairs of regions were included if they were at least 20mm apart in Euclidean distance. (c-d) Control analysis considering region pairs with at least 50 stimulation experiments (probability) or significant responses (amplitude). (e-f ) Control analysis in which pairs of regions were considered to be anatomically connected (unconnected) if they shared a connection in at least 90% (no more than 10%) of participants. (g-h) Control analysis in which response to stimuli were considered significant for z-scored CCEP amplitude larger than 4. (i-j) Control analysis in which response to stimuli were considered significant for z-scored CCEP amplitude larger than 6. (k-l) Control analysis in which no restrictions were imposed on Euclidean distance between pairs of regions.

- 68 ;

 68 Fig 1e) and amplitude (p < 10 -12 ; Fig 1i), compared to unconnected ones. As shown in Figs 1f

FIG. 1 .

 1 FIG. 1. (a)Schematics of the F-TRACT project. Subjects were implanted with SEEG electrodes capable of delivering localized stimuli and recording electrical activity at millisecond resolution. The toy example shows a stimulation performed on a single subject, where the effects of a single-pulse stimulus are recorded (schematics are illustrative and do not reflect the bipolar stimulation montage and the true spacing between contacts). A significant response is measured by contact 1, as defined by a z-scored CCEP larger than 5. Contact 2 does not record a significant response. Based on its spatial coordinate, each contact was assigned to a cortical gray matter region. Whole-brain matrices of response probability (P) and median amplitude (A) were computed by combining data from 29,055 stimulations in 550 subjects. (b) Normative, group-representative measures of brain connectivity were inferred from 1,000 healthy participants of the Human Connectome Project using structural, diffusion, and resting-state functional MRI data. (c) Pairs of regions used to model response probability. Violet and copper matrix entries denote regions that are likely and unlikely to be anatomically connected, respectively. Warm-and cool-colored line plots indicate the number of times each region featured as a stimulation and recording site, respectively, and were projected onto the cortical surface on panel (d). (e) Response probability distributions for all (black), connected (violet), and unconnected (copper) region pairs. (f ) Bars show the absolute value of Spearman's ρ between response probability and ED, SC and FC, computed for all (left), connected (center), and unconnected (right) region pairs. Circles and error bars indicate the mean and 95% confidence interval correlations obtained from 1,000 bootstrapped samples. (g-j) Same as (c-f) but for response amplitude.

FIG. 2 .

 2 FIG. 2. (a) A significant response is observed between anatomically unconnected stimulation and recording sites. How does white matter connectivity facilitate and constrain stimulus propagation? (b) Schematics of five putative measures of brain network communication. Briefly, shortest path efficiency (SPE) quantifies communication via the least costly path between regions; navigation efficiency (NE) presumes propagation via geometrically greedy paths; search information (SI) assumes a diffusive communication policy and measures the opportunity for signal dispersion along shortest paths; communicability (CMY) is a diffusive broadcasting process; and diffusion efficiency (DE) models propagation via random walks. (c) Network communication measures rely on different aspects of structural connectivity to quantify inter-areal propagation. Left: SPE models stimulus response (R) between stimulation (A, B, C ) and recording (T ) sites based exclusively on properties of the shortest-e.g., fewest intermediate regions-path between them. Center: SI takes into account the local connectivity along shortest paths. A stimulus propagated via regions with many connections is more prone to dispersion, leading to weaker responses. Despite being connected to T via a single intermediate region, A, B, and C, in this order, provide decreasing opportunities for signal dispersion. Right: Communicability is a broadcasting process that unfolds through the entire network. Short paths contribute more to communication than long ones. Despite both being connected to T via a single intermediate region, stimulation to A evokes stronger responses than stimulation to B, since A has a larger number of short paths to T.

  FIG. 3. (a-b) Spearman's correlations to response probability for (a) all, (b) anatomically connected and (c) anatomically unconnected region pairs. Absolute values of full and partial (controlling for ED) correlations are shown in light-and darkcolored bars, respectively. Gray bars indicate correlations that were not significant at α = 5%. Circles and error bars indicate, respectively, the mean and 95% confidence interval correlations obtained from 1,000 bootstrapped samples. Black triangles mark full correlations illustrated by scatter plots in the lower panels. Each data point in the scatter plots represents a pair of regions, with darker colors indicating higher point density. (d-f ) Same as (a-b) but for response amplitude.

FIG. 4 .

 4 FIG. 4. Machine learning predictions of stimulus propagation. A lasso regression model was trained to predict stimulus response from a combination of SC, FC, ED, and all network communication measures. By penalizing model complexity (number of predictors), the lasso regression seeks to identify a sparse set of input variables that accurately predicts the outcome variable. (a) Out-of-sample prediction accuracy for 15 repetitions of 10-fold cross validation, computed as the Pearson correlation between predicted and observed probabilities. (b) Scatter plot of predicted versus observed probabilities with prediction accuracy evaluated as the Pearson correlation, R 2 , and root mean squared error (RMSE). (c) Average contribution of input predictors across 15 repetitions of 10-fold cross validation. (d-f ) Same as (a-c) but for predictions of response amplitude.

  is supported by a Melbourne Research Scholarship. AZ is supported by the National Health and Medical Research Council of Australia (APP1118153).

  FIG. S1. Pairwise t-statistics from two-sample t-tests performed on the distributions of 1,000 bootstrapped Spearman correlation coefficients presented in Fig 3. A warm-colored (cool-colored) matrix entry ij indicates that the mean correlation obtained for measure i is statistically higher (lower) than the mean correlation obtained for measure j. Matrices are symmetric about their main diagonal, with the opposite signs. A white matrix entry indicates no significant difference between the mean correlations obtained from two measures. Statistical tests for each n × n t-statistics matrix were adjusted using the Benjamini-Hochberg correction of the False Discovery Rate for (n × (n -1))/2 multiple comparisons. Statistics are shown for full Spearman correlations with response probability for (a) all, (b) anatomically connected, and (c) anatomically unconnected region pairs. (d-e) Same as (a-c) but for response amplitude. (g-l) Same as (a-f) but for partial Spearman correlations controlling for the effect of Euclidean distance.

ACKNOWLEDGMENTS

Brain imaging data were provided by the Human Connectome Project, WU-Minn Consortium (1U54MH091657; Principal Investigators: David Van Essen and Kamil Ugurbil) funded by the 16 National Institutes of Health (NIH) institutes and centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. CS was supported by the Australian Research Council (grant number DP170101815). MJ and OD were supported by funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement no. 616268 F-TRACT, the European Union's Horizon 2020 Framework Programme for Research and Innovation under Specific Grant Agreement No. 785907 and 945539 (Human Brain Project SGA2 and SGA3), and from the Agence Nationale de la Recherche grant number ANR-21-NEUC-0005-01. SML


R Ravaud

Senior IEEE G Lemarquand
email: guy.lemarquand@ieee.org

V Lemarquand

Force and Stiffness of Passive Magnetic Bearings Using Permanent Magnets. Part 2 : Radial Magnetization

Keywords: Magnetic forces, analytical calculation, ring permanent magnet, magnetic bearing I

This paper deals with the calculation of the force and the stiffness between two ring permanent magnets whose polarization is radial. Such a configuration corresponds to a passive magnetic bearing.

The magnetic force exerted between ring permanent magnets is determined by using the coulombian model. The expressions obtained are semi-analytical and we show that it is not possible to find an exact analytical expression of the force between two ring permanent magnets. Then, thanks to these semianalytical calculations, the ring dimensions are optimized in order to have a great force or a great stiffness.

Moreover, we show that the relative position of the rings for which the force is the strongest depends on the air gap dimension. This result is new because the curvature effect is taken into account in this paper.

We can say that such semi-analytical expressions are more precise than the numerical evaluation of the magnetic forces obtained with the finite element method. Moreover, semi-analytical expressions have a low computational cost whereas the finite element method has a high one. Thereby, as shown in this paper, such calculations allow an easy optimization of quadripolar lenses or devices using permanent magnets.

Permanent magnets are commonly used in many electrical devices and engineering applications. Most engineering applications need several ring permanent magnets and the determination of the magnetic force between them is thus required. It is to be noted that several way of obtaining the magnetic forces between rings are possible. Authors generally use either numerical means or 2D analytical calculations for determining magnetic fields or magnetic forces created between ring permanent magnets [START_REF] Azzerboni | Computation of the magnetic field in massive conductor systems[END_REF]- [START_REF] Akoun | 3d analytical calculation of the forces exerted between two cuboidal magnets[END_REF]. One of the major problem of a numerical approach lies in the fact that it has a high computational cost.

Moreover, numerical algorithms are not so precise as analytical calculations. Alternative solutions are thus required. As it is difficult to obtain fully analytical expressions of the magnetic field created by ring permanent magnets [START_REF] Kim | Mutual inductance of noncoaxial circular coils with constant current density[END_REF]- [START_REF] Kim | Restoring force between two noncoaxial circular coils[END_REF], the determination of the forces between them is still more difficult.

Consequently, authors use semi-analytical expressions which are very important steps in order to evaluate either the magnetic field created by ring permanent magnets or the magnetic forces exerted between them.

The determination of the magnetic field created by ring permanent magnets have been studied a lot with such approaches [START_REF] Babic | New expressions for calculating the magnetic field created by radial current in massive disks[END_REF]- [START_REF] Rakotoarison | Using coulombian approach for modeling scalar potential and magnetic field of a permanent magnet with radial polarization[END_REF].

Another way of calculating the forces between ring permanent magnets can be done by using the 2D analytical expressions of the magnetic field created by infinite parallelepipedic magnets. Assuming that the magnetic field created by ring permanent magnets can be approximated by the magnetic field created by infinite parallelepipedic magnets either to the faces of the magnets (2D analytical approach) [START_REF] Lemarquand | Nonsinusoidal torque of permanent-magnet couplings[END_REF]- [START_REF] Zhu | Comparison of halbach magnetized brushless machines based on discrete magnet segments or a single ring magnet[END_REF], the magnetic forces can be calculated more easily with the 2D analytical approach.

However, these formula are not valid when the ring radius is small [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF]- [START_REF] Ravaud | Magnetic field produced by a tile permanent magnet whose polarization is both uniform and tangential[END_REF]. In any case, the magnetic field created by ring permanent magnets can be determined in terms of elliptic integrals [START_REF] Durand | Magnetostatique[END_REF] for ring permanent magnets axially magnetized. Such an approach is appropriate because the algorithms used to calculate elliptic integrals are both very robust and fast. Consequently, the expressions obtained have a very low computational cost (less than 0.2 s to determine the magnetic components of the field created by ring permanent magnets whose polarization is axial).

First, this paper presents useful semi-analytical expressions of the force and the stiffness exerted between two ring permanent magnets whose polarization is radial. Then, this paper explains why we cannot reduce some numerical integrations of the semi-analytical expressions of the force between two ring permanent magnets. Eventually, this paper deals with the optimization of the ring dimensions in order to have either a great force or a great stiffness between the rings. We show that the relative position between the rings is of great importance in the design of devices using ring permanent magnets. All the expressions determined 1. Representation of the two ring permanent magnets which are radially centered. The inner radius of the outer ring is denoted r in , the outer one is rout, its height is h. The inner radius of the inner ring is denoted r in2 , the outer one is r out2 , its height is

z b -za.
in this paper are available online [30].

II. CALCULATION OF THE AXIAL FORCE BETWEEN TWO RING PERMANENT MAGNETS WHOSE

POLARIZATION IS RADIAL

This section presents a semi-analytical calculation of the force exerted between two ring permanent magnets whose polarization is radial and which are radially centered. Such configuration corresponds to an axial passive magnetic bearing using ring magnets radially magnetized for mutual attraction. We can say that the devices realized with ring permanent magnets radially magnetized were the first to be built. Moreover, they are usually made of several sections which are stacked together.

A. Notation and geometry

The geometry considered is shown in Fig 1 . A two dimensional representation of the passive magnetic bearing is shown in Fig 2.

The outer radius of the outer ring is r out and the inner one is r in . The outer ring height is h. The outer radius of the inner ring is r out2 and the inner one is r in2 . The inner ring height is z bz a . It is to be noted . Radial magnetic bearing using magnets radially magnetized for mutual attraction. The inner radius of the outer ring is denoted r in , the outer one is rout, its height is h. The inner radius of the inner ring is denoted r in2 , the outer one is r out2 , its

height is z b -za.
that the coulombian model of a permanent magnet is used. Consequently, each ring permanent magnet is represented by both two curved planes which correspond to the inner and outer faces of the rings which are charged with a surface magnetic pole density σ * and a magnetic pole volume density σ * v . For each case, the inner face is charged with the surface magnetic pole density +σ * and the outer one is charged with the surface magnetic pole density -σ * . It is noted that all the illustrative calculations are done with σ * = J. n = 1T where J is the magnetic polarization vector and n is the unit normal vector which is directed towards 0. Moreover, it is noted that the magnetic pole volume density exits for ring permanent magnets whose polarization is radial in order to to have a charge equilibrium of the ring magnet.

B. Semi-analytical expression of the magnetic force

The two ring permanent magnets which form an axial passive bearing are supposed to be radially centered. Consequently, there is only the axial component of the magnetic force which is exerted between the two rings. We can call it the axial force F z . This axial force can be determined by integrating the magnetic field created by the outer ring on the contributions of the inner one. We must take into account both the magnetic pole surface density and the magnetic pole volume density of each ring. Consequently, as there are two magnetic pole surface densities and one magnetic pole volume density for each ring permanent magnet, we have nine terms for determining the axial force between the two ring permanent magnets. By denoting H z the axial component of the magnetic field produced by the outer ring permanent magnet, Thus, this axial force F z can be written as follows:

F z = (Sin) H z σ * 2 d S - (Sout) H z σ * 2 d S + (V ) H z σ * 2 r 2 d Ṽ (1)
where σ * 2 is the magnetic pole surface density owing to the inner ring permanent magnet, (S in ) is the inner face of the inner ring permanent magnet and (S out ) is the outer face of the inner ring permanent magnet and (V ) is the volume of the inner ring permanent magnet. Therefore, the axial force F z can be expressed as follows:

F z = - 2π θ1=0 h z1=0 2π θ2=0 z b z2=za a(r in , r out2 )dz 1 dθ 1 dz 2 dθ 2 + 2π θ1=0 h z1=0 2π θ2=0 z b z2=za a(r out , r out2 )dz 1 dθ 1 dz 2 dθ 2 + 2π θ1=0 h z1=0 2π θ2=0 z b z2=za a(r in , r in2 )dz 1 dθ 1 dz 2 dθ 2 - 2π θ1=0 h z1=0 2π θ2=0 z b z2=za a(r out , r in2 )dz 1 dθ 1 dz 2 dθ 2 - rout r1=rin 2π θ1=0 h z1=0 2π θ2=0 z b z2=za b(r out2 , r 1 )dr 1 dz 1 dθ 1 dz 2 dθ 2 + rout r1=rin 2π θ1=0 h z1=0 2π θ2=0 z b z2=za b(r in2 , r 1 )dr 1 dz 1 dθ 1 dz 2 dθ 2 + rout2 r2=rin2 2π θ1=0 h z1=0 2π θ2=0 z b z2=za b(r in , r 2 )dz 1 dθ 1 dr 2 dz 2 dθ 2 - rout2 r2=rin2 2π θ1=0 h z1=0 2π θ2=0 z b z2=za b(r out , r 2 )dz 1 dθ 1 dr 2 dz 2 dθ 2 + rout r1=rin rout2 r2=rin2 2π θ1=0 h z1=0 2π θ2=0 z b z2=za c(r 1 , r 2 )dr 1 dz 1 dθ 1 dr 2 dz 2 dθ 2 (2) with a(α, β) = σ * 1 σ * 2 4πµ 0 (z 2 -z 1 )αβ (α 2 + β 2 -2αβ + (z 2 -z 1 ) 2 ) 3 2 (3)
where σ * 1 is the magnetic pole surface density owing to the inner ring permanent magnet. and

b(α, β) = a(α, β) β (4) c(α, β) = a(α, β) αβ (5)
The next step is to evaluate (2). If we determined directly F z for given values with all the integrals presented in (2), the computational would be too high because several numerical integrals should be determined at a time. Consequently, we must reduce the number of integrals by integrating analytically a(α, β), b(α, β) and c(α, β) according to the integral variables. Unfortunately, it is not possible to find a fully analytical expression of the magnetic force between two ring permanent magnets but we can use a useful semi-analytical expression which is expressed as follows:

F z = σ * 1 σ * 2 2µ 0 2π θ1=0 Sdθ 1 + σ * 1 σ * 2 2µ 0 2π θ1=0 M dθ 1 + σ * 1 σ * 2 2µ 0 2π θ1=0 rout2 r2=rin2 V dθ 1 dr 2 (6)
where S corresponds to the magnetic interaction between the magnetic pole surface densities of each ring permanent magnet. Then, M corresponds to the magnetic interaction between the magnetic pole surface densities of one ring permanent magnet and the magnetic pole volume density of the other one.

At least, V corresponds to the magnetic interaction between the magnetic pole volume densities of each ring permanent magnet.

The first contribution S is given by [START_REF] Selvaggi | Calculating the external magnetic field from permanent magnets in permanent-magnet motors -an alternative method[END_REF].

S = -f (z a , z b , h, θ 1 , r in , r out2 ) +f (z a , z b , h, θ 1 , r out , r out2 ) +f (z a , z b , h, θ 1 , r in , r in2 ) -f (z a , z b , h, θ 1 , r out , r in2 ) (7)
with

f (α 1 , α 2 , α 3 , θ 1 , α 5 , α 6 ) = α 5 α 6 log α 3 -α 1 + α 2 5 + α 2 6 + (α 3 -α 1 ) 2 -2α 5 α 6 cos(θ 1 ) +α 5 α 6 log α 1 + α 2 5 + α 2 6 + α 2 1 -2α 5 α 6 cos(θ 1 ) -α 5 α 6 log α 3 -α 2 + α 2 5 + α 2 6 + (α 3 -α 2 ) 2 -2α 5 α 6 cos(θ 1 ) -α 5 α 6 log α 2 + α 2 5 + α 2 6 + α 2 2 -2α 5 α 6 cos(θ 1 ) (8)
The second contribution M is given by ( 9).

M = -t(r in , r out , r out2 , h, z a , z b , θ 1 ) +t(r in , r out , r in2 , h, z a , z b , θ 1 ) +t(r in , r out , r in , h, z a , z b , θ 1 ) -t(r in , r out , r out , h, z a , z b , θ 1 ) (9)
with

t(β 1 , β 2 , β 3 , β 4 , β 5 , β 6 , θ 1 ) = t (2) β 1 , β 2 , β 4 -β 5 , β 2 3 + (β 4 -β 5 ) 2 , 2β 3 cos(θ 1 )
+t (2)

β 1 , β 2 , β 5 , β 2 3 + β 2 5 , 2β 3 cos(θ 1 ) -t (2) β 1 , β 2 , β 4 -β 6 , β 2 3 + (β 4 -β 6 ) 2 , 2β 3 cos(θ 1 ) -t (2) β 1 , β 2 , β 6 , β 2 3 + β 2 6 , 2β 3 cos(θ 1 ) (10)
and

t (2) (β 1 , β 2 , q, d, f ) = t (3) (β 2 , q, d, f ) -t (3) (β 1 , q, d, f ) (11)
and

t (3) (s, q, d, f ) = -s + 4d -f 2 -4q 2 2 arctan -f + 2s 4d -f 2 -4q 2 - f 4 log d -q 2 -f s + s 2 +s log q + d -f s + s 2 + q log -f + 2(s + d -f s + s 2 ) - 4d -f 2 -4q 2 + f η log [u 1 ] 4η - -4d + f 2 + 4q 2 + f η log [u 2 ] 4η (12)
u 1 = - 2 f 2 + 4f q 2 -f 2 (η + 2s) + 4d(-f + η + 2s) q 2 (-4d + f 2 + 4q 2 -f η)(-f + η + 2s) - 8q(-2qs + η d -f s + s 2 ) q 2 (-4d + f 2 + 4q 2 -f η)(-f + η + 2s) (13)
u 2 = - 2 f 2 + 4f q 2 -f 2 (η -2s) -4d(f + η -2s) q 2 (-4d + f 2 + 4q 2 + f η)(f + η -2s) - -8q(2qs + η d -f s + s 2 ) q 2 (-4d + f 2 + 4q 2 + f η)(f + η -2s) (14)
with

η = -4d + f 2 + 4q 2 (15)
The third contribution V is given by ( 16).

V = th (1) (r out , r 2 , z a , z b , h, θ 1 ) -th (1) (r in , r 2 , z a , z b , h, θ 1 ) (16)
with

th (1) = t (3) (r 1 , h -z a , r 2 2 + (h -z a ) 2 , 2r 2 cos(θ 1 )) +t (3) (r 1 , z a , r 2 2 + z 2 a , 2r 2 cos(θ 1 )) -t (3) (r 1 , h -z b , r 2 2 + (h -z b ) 2 , 2r 2 cos(θ 1 )) -t (3) (r 1 , z b , r 2 2 + z 2 b , 2r 2 cos(θ 1 )) (17)

C. Expression of the axial stiffness between two rings whose polarization is radial

The stiffness K exerted between two ring permanent magnets can be determined by calculating the derivative of the axial force with respect to z a . We set z b = z a + b where b is the height of the inner ring permanent magnet. Thus, the axial stiffness K can be calculated with [START_REF] Blache | New structures for linear displacement sensor with hight magnetic field gradient[END_REF].

K = - ∂ ∂z a F z ( 18
)
where F z is given by (2). We obtain :

K = K S + K M + K V ( 19
)
where K S corresponds to the stiffness determined with only the magnetic pole surface densities of each ring permanent magnet. Then, K M corresponds to the stiffness determined with the magnetic interaction between the magnetic pole surface densities of one ring permanent magnet and the magnetic pole volume density of the other one. At least, K V corresponds to the stiffness determined with the magnetic interaction between the magnetic pole volume densities of each ring permanent magnet. Thus, the first contribution K S is expressed as follows:

K S = η 31 1 √ α 31 K * - 4r 3 r 1 α 31 - 1 √ β 31 K * - 4r 3 r 1 β 31 + 1 √ δ 31 K * - 4r 3 r 1 δ 31 - 1 √ γ 31 K * - 4r 3 r 1 γ 31 +η 41 1 √ α 41 K * - 4r 4 r 1 α 41 - 1 √ β 41 K * - 4r 4 r 1 β 41 + 1 √ δ 41 K * - 4r 4 r 1 δ 41 - 1 √ γ 41 K * - 4r 4 r 1 γ 41 +η 32 1 √ α 32 K * - 4r 3 r 2 α 32 - 1 √ β 32 K * - 4r 3 r 2 β 32 + 1 √ δ 32 K * - 4r 3 r 2 δ 32 - 1 √ γ 32 K * - 4r 3 r 2 γ 32 +η 42 1 √ α 42 K * - 4r 4 r 2 α 42 - 1 √ β 42 K * - 4r 4 r 2 β 42 + 1 √ δ 42 K * - 4r 4 r 2 δ 42 - 1 √ γ 42 K * - 4r 4 r 2 γ 42 (20)
with

η ij = 2r i r j σ * µ 0 (21)
α ij = (r i -r j ) 2 + z 2 a ( 22
)
β ij = (r i -r j ) 2 + (z a + h) 2 (23)
γ ij = (r i -r j ) 2 + (z a -h) (24)
δ ij = (r i -r j ) 2 + (b -h) 2 + z a (2b -2h + z a ) (25)
K * [m] = π 2 0 1 1 -m sin(θ) 2 dθ ( 26
)
The second contribution K M is expressed as follows:

K M = σ * 1 σ * 2 2µ 0 2π θ=0 udθ (27)
with

u = f (r in , r out , r in2 , h, z a , b, θ) -f (r in , r out , r out2 , h, z a , b, θ) +f (r in2 , r out2 , r in , h, z a , b, θ) -f (r in2 , r out2 , r out , h, z a , b, θ) (28)
and

f (α, β, γ, h, z a , b, θ) = -γ log α -γ cos(θ) + α 2 + γ 2 + z 2 a -2αγ cos(θ) +γ log α -γ cos(θ) + α 2 + γ 2 + (z a + b) 2 -2αγ cos(θ) +γ log α -γ cos(θ) + α 2 + γ 2 + (z a -h) 2 -2αγ cos(θ) -γ log α -γ cos(θ) + α 2 + γ 2 + (b -h) 2 + 2z a (b -h) + z 2 a -2αγ cos(θ) +γ log β -γ cos(θ) + β 2 + γ 2 + z 2 a -2αγ cos(θ) -γ log β -γ cos(θ) + β 2 + γ 2 + (z a + b) 2 -2αγ cos(θ) +γ log β -γ cos(θ) + β 2 + γ 2 + (z a -h) 2 -2αγ cos(θ) -γ log β -γ cos(θ) + β 2 + γ 2 + (b -h) 2 + 2z a (b -h) + z 2 a -2αγ cos(θ) (29)
The third contribution K V is expressed as follows:

K V = σ * 1 σ * 2 2µ 0 2π θ=0 rout r1=rin δdθ (30)
with

δ = -log r in2 -r 1 cos(θ) + r 2 1 + r 2 in2 + z 2 a -2r 1 r in2 cos(θ) + log r in2 -r 1 cos(θ) + r 2 1 + r 2 in2 + (z a + b) 2 -2r 1 r in2 cos(θ) -log r in2 -r 1 cos(θ) + r 2 1 + r 2 in2 + (b -h) 2 + 2bz a -2hz a + z 2 a -2r 1 r in2 cos(θ) + log r in2 -r 1 cos(θ) + r 2 1 + r 2 in2 + (z a -h) 2 -2r 1 r in2 cos(θ) + log r out2 -r 1 cos(θ) + r 2 1 + r 2 out2 + z 2 a -2r 1 r out2 cos(θ) -log r out2 -r 1 cos(θ) + r 2 1 + r 2 out2 + (z a + b) 2 -2r 1 r out2 cos(θ) -log r out2 -r 1 cos(θ) + r 2 1 + r 2 out2 + (z a -h) 2 -2r 1 r out2 cos(θ) + log r out2 -r 1 cos(θ) + r 2 1 + r 2 out2 + (b -h) 2 + 2bz a -2hz a + z 2 a -2r 1 r out2 cos(θ) (31)
We can say that the expression of the axial stiffness can be determined analytically if we only take into account the magnetic pole surface densities of each ring.

III. DISCUSSION ABOUT THE POSSIBILITY OF REDUCING THE NUMBER OF NUMERICAL INTEGRATIONS FOR THE AXIAL FORCE EXPRESSION

The aim of this section is to explain why we cannot reduce the number of numerical integrations of the semi-analytical expression of the force exerted between two ring permanent magnets radially magnetized.

They correspond physically to three kinds of interactions between two ring permanent magnets whose polarization is radial.

A. Interaction between the contributions of the surface densities

The first kind of physical interaction is due to the surface contributions of each ring permanent magnet.

This physical interaction corresponds to the case when the surface densities of the outer ring are integrated with the surface densities of the inner one. In Eq. ( 2), these surface contributions correspond to the integrand denoted a(α, β) where α and β can be r in , r out , r in2 and r out2 . Let us consider the integrand a(α, β). As we can see in Eq. ( 2), the integration variables depend on θ 1 , θ 2 , z 1 and z 2 . The integration according θ 2 does not change the form of the integrand a(α, β) because a(α, β) does not depend on θ 2 .

Consequently, we can say that the form of the force F s between the surface density contributions of each ring permanent magnet is given as follows :

F s = θ1 z1 z2 a 1 (z 2 -z 1 ) (a 2 -a 3 cos(θ 1 ) + (z 2 -z 1 ) 2 ) 3 2 dθ 1 dz 1 dz 2 (32)
where a 1 , a 2 and a 3 are constant. For example, we can have

a 1 = 2π(h -z a ), a 2 = 2π(r 2 in + r 2 out + (h -z a )
2 ) and a 3 = 2π(-2r in r out2 ) for our illustration here. In short, these parameters are given by [START_REF] Yonnet | Permanent magnet bearings and couplings[END_REF] and correspond to the case when four integrals must be determined.

where α and β can be r 1 or r 2 . Let us consider the integrand c(α, β). As we can see in Eq. ( 2), the integration variables depend on θ 1 , θ 2 , z 1 , z 2 , r 1 and r 2 . The integration according θ 2 does not change the form of the integrand c(α, β) because c(α, β) does not depend on θ 2 . Consequenlty, we can say that the form of the force F v between the volume density contributions is given as follows :

F v = θ1 r1 r2 z1 z2 (z 2 -z 1 ) (c 1 -c 2 cos(θ 1 ) + (z 2 -z 1 ) 2 ) 3 2 dθ 1 dz 1 dz 2 dr 1 dr 2 (35)
where c 1 and c 2 are constant. These parameters are given by ( 2) and correspond to the case when six integrals must be determined.

After having integrated according z 1 , z 2 and r 1 , we obtain the same integrand as the one obtained previously in the case of the study of the force exerted between the surface and volume charge contributions.

Consequently, we deduct that we cannot integrate analytically F v according to θ 1 . Moreover, the analytical integration according to r 2 does not seem possible. As a conclusion, we can say that the obtaining of a fully analytical expression of the force between two ring permanent magnets whose polarization is radial does not seem possible but a semi-analytical expression can be used to determine this axial force.

IV. OPTIMIZATION OF THE INNER RING PERMANENT MAGNET DIMENSIONS

This section discusses the optimal dimensions of the rings in order to have either a great force or a great stiffness.

A. Influence of the air gap dimension on the force and the stiffness

First, we study the influence of the air gap dimension on the force and the stiffness between the two ring permanent magnets. For this purpose, we represent the axial force versus the axial displacement of the inner ring for different air gaps in Fig. 3. It is noted that, in our configuration, the air gap corresponds to the difference between r in and r out2 . Furthermore, the width and the height of each ring permanent magnet are constant. We take r in = 0.025m, r out = 0.028m, r out2 -r int2 = 0.003m, J = 1T, h = 0.003m, zbz a = 0.003m.

Fig. 3 shows three important points.

First, we see that the smaller the air gap between the ring permanent magnets is, the greater the axial force is. Consequently, it is necessary to have the smaller air gap between ring permanent magnets if a great force is searched. This result is well-known. It was shown with the two dimensional approach.

Second, we see that the exact position of the maximal force exerted between two ring permanent magnets Fig. 3. Representation of the axial component of the magnetic force exerted between two ring permanent magnets versus the axial displacement of the inner ring permanent magnet for different air gaps; r in = 0.025m, rout = 0.028m,r out2 -r int2 = 0.003m,

J = 1T, h = 0.003m, zb -za = 0.003m.
depends slightly on the air gap dimension. This result is new because our study uses a three-dimensional approach of the magnetic force whereas the previous ones used a two-dimensional approach. The magnet curvature must be taken into account in order to obtain precisely the position of the maximal force exerted between two ring magnets. The exact position of the maximal force is represented in

B. Determination of the optimal height of the inner ring permanent magnet

Another parameter which can be optimized is the height of the inner ring permanent magnet. To do so, the axial component of the magnetic force is plotted versus the axial displacement of the inner ring Sti ness [N/m] Fig. 5. Representation of the axial stiffness of the magnetic force exerted between two ring permanent magnets versus the axial displacement of the inner ring permanent magnet for different air gaps; r in = 0.025m, rout = 0.028m,r out2 -r int2 = 0.003m, J = 1T, h = 0.003m, zbza = 0.003m. permanent magnet for several inner ring heights in Fig. 6. The values taken for the parameters are still the same as the previous ones. Fig. 6 shows that the axial component of the magnetic force is the greatest if its height equals the outer ring height. Indeed, if the ring inner height is smaller than the ring outer one, the smaller the inner ring width is, the smaller the axial component of the magnetic force is. If the ring inner height is greater than the ring outer one, the greater the inner ring width is, the smaller the axial component of the magnetic force is. Consequently, if a great force is searched, the two ring permanent magnets must have the same height. Furthermore, we see that the gradient of the curves is the more important when the inner ring height is the same as the outer one. This result is consistent with Fig 7 where the axial stiffness is represented versus the axial displacement of the inner ring permanent magnet.

C. Determination of the optimal width of the inner ring permanent magnet

The third parameter which can be optimized in our configuration is the width of the inner ring permanent magnet. To do so, the axial component of the magnetic force is plotted versus the axial displacement of the inner ring permanent magnet for several inner ring widths in Fig. 8. The values taken for the parameters are the same as the previous ones. Fig. 8 shows that the greater the inner ring width is, the greater the axial component of the magnetic force is. However, it is noted that a compromise in the ring dimensions must be found because the cost of the magnet must be taken into account. A good compromise can be found as follows : if the inner ring width equals two times its height, the axial component of the magnetic force is 72N when z = 0.0015m. If the inner ring width equals three times its height, the axial component of the magnetic force is 73N. Consequenlty, we deduct that it is not necessary to have an inner ring width which is greater than two times its height.

The optimal stiffness depends also on the inner ring width. To see that, we have represented in Fig. 9. Representation of the axial stiffness exerted between two ring permanent magnets versus inner ring width r in = 0.025m, rout = 0.028m,r out2 = 0.0249m, J = 1T, h = 0.003m, zbza = 0.003m.

is, the greater the axial stiffness is.

V. OBTAINING THE BEST CONFIGURATION

The previous section shows that the ring dimensions must be optimized in order to create a very good passive magnetic bearing. The air gap must be the smallest, the ring heights must be the same and we have shown that the inner ring width must equal two times its height. All these parameters have been determined with a outer ring whose cross-section is a square. However, we can also optimize the outer ring in order to improve the passive bearing. By taking into account the optimal dimensions found in the previous section, we can compare three configurations. The first one, shown in Consequently, this last parameter must be taken into account in the design of passive bearings using ring permanent magnets.

C A B Fig. 10. Representation of three passive magnetic bearings. A = r in = 0.025m, rout = 0.028m, r out2 = 0.0249m, r in2 = 0.0219m, h = 0.003m, zbza = 0.003m. B =r in = 0.025m, rout = 0.031m, r out2 = 0.0249m, r in2 = 0.0189m, h = 0.003m, zbza = 0.003m. C= r in = 0.025m, rout = 0.028m, r out2 = 0.0249m, r in2 = 0.0219m, h = 0.006m, zbza = 0.006m.

VI. CONCLUSION

This paper has presented new three-dimensional semi-analytical expressions allowing us to determine both the axial force and the axial stiffness between two ring permanent magnets whose polarization is radial. This paper also discusses the reason why we cannot find a fully analytical expression of the axial component of the magnetic force. Eventually, we have discussed the optimal dimensions of the ring permanent magnets which allow us to have either a great axial force or a great axial stiffness. We have shown that a good compromise can be found when the cross-section of a ring is a rectangle. However, the relative position of the two rings is not constant but depends on the air gap dimension. Such results can be very useful for people involved in the design of magnetic bearings.

  Fig. 1. Representation of the two ring permanent magnets which are radially centered. The inner radius of the outer ring is denoted

  Fig.2. Radial magnetic bearing using magnets radially magnetized for mutual attraction. The inner radius of the outer ring is

Fig 4 .

 4 Such result is very useful because it clearly shows that if a great axial force is searched, the relative height between the two rings inner ring depends on the air gap dimension.Eventually, Fig3shows that the stiffness depends greatly on the air gap dimension. Fig3shows that the smaller the air gap dimension is, the greater the stiffness is because the gradient of the curve is the more important for small air gaps. This result is consistent with the representation of the axial stiffness versus the axial displacement of the inner ring permanent magnet (Fig5). Indeed, we see that the smaller the air gap dimension is, the greater the axial stiffness is. Moreover, we see that when the axial force is maximal in Fig 3, the axial stiffness equals zero in Fig 5, which is still consistent.

Fig. 4 .

 4 Fig.[START_REF] Furlani | A three-dimensonal field solution for radially polarized cylinders[END_REF]. Representation of the position of the maximal value of the axial force versus the air gap dimension; r in = 0.025m, rout = 0.028m,r out2 -r int2 = 0.003m, J = 1T, h = 0.003m, zbza = 0.003m.

Fig. 6 .

 6 Fig.6. Representation of the axial component of the magnetic force exerted between two ring permanent magnets versus the axial displacement of the inner ring permanent magnet for different heights; r in = 0.025m, rout = 0.028m,r out2 = 0.0249m, J = 1T, h = 0.003m, za = 0m.

Fig. 7 .Fig. 8 .

 78 Fig.7. Representation of the axial component of the magnetic force exerted between two ring permanent magnets versus the axial displacement of the inner ring permanent magnet for different heights; r in = 0.025m, rout = 0.028m,r out2 = 0.0249m, J = 1T, h = 0.003m, za = 0m.

  Fig 10-A, consists of two rings whose cross-section is a square. The second one, shown in Fig 10-B, consists of two rings whose cross-section is a rectangle whose width equals two times its height . The third one, shown in Fig 10-C, consists of two rings whose cross-section is a rectangle whose height equals two times its width. For each structure presented in Fig 10, the axial force and the axial stiffness are determined versus the axial displacement of the inner ring permanent magnet. The axial force is shown in Fig 11 and the axial stiffness is shown in Fig 12.

Figs 11 and 12

 12 Figs 11 and 12 show that the best configration is the configuration B presented in Fig 10. However, the relative height between the two ring permanent magnets depends on the air gap dimension (Fig 4).

Fig. 11 .Fig. 12 .

 1112 Fig.11. Representation of the axial force exerted between two ring permanent magnets versus inner ring width r in = 0.025m, rout = 0.028m,r out2 = 0.0249m, J = 1T, h = 0.003m, zbza = 0.003m.

After having integrated according z 1 and z 2 , we obtain a semi-analytical expression with only one numerical integration whose variable depends on θ 1 :

These parameters are given by ( 8) and they are not a function of the angle θ. It is to be noted that (33) cannot be integrated analytically.

B. Interaction between the contributions of the surface and volume densities

The second kind of physical interaction is due to both the surface and the volume contributions of each ring permanent magnet. This physical interaction corresponds to the case when the surface densities of each ring permanent are integrated with the volume densities of the other ones. In Eq. ( 2), these surface contributions correspond to the integrand denoted b(α, β) where α and β can be r in , r out , r in2 , r out2 , r 1 and r 2 . Let us consider the integrand b(α, β). As we can see in Eq. ( 2), the integration variables depend on θ 1 , θ 2 , z 1 , z 2 , r 1 and r 2 . The integration according θ 2 does not change the form of the integrand b(α, β) because b(α, β) does not depend on θ 2 . Consequently, we can say that the form of the force F vs between the surface density contributions and the volume density contributions of each ring permanent magnet is given as follows :

where b 1 , b 2 and b 3 are constant. These parameters are given by ( 2) and correspond to the case when five integrals must be determined.

After having integrated according z 1 , z 2 and r 1 , we obtain a semi-analytical expression with only one numerical integration whose variable depends on θ 1 : this semi-analytical expression is in fact given by (33) where s depends on θ 1 . We see that (34) owns a term which has the same form as the one presented in Eq.( 33). Consequently, it cannot be integrated analytically as well.

C. Interaction between the contributions of the volume densities

The third kind of physical interaction is due to the volume contributions of each ring permanent magnet.

This physical interaction corresponds to the case when the volume densities of each ring permanent magnet are integrated together. In Eq. ( 2), this volume contribution correspond to the integrand denoted c(α, β)