Jordan Orly

MD Annouk Bisdorff

MD Antoine Fraissenon

MD Aline Joly

MD Grégoire Boulouis

PhD Laurent Guibaud

MD, PhD; Elsa Tavernier

PhD Stéphanie Mallet

Clément Marcelin

MD; Juliette Miquel

Ludovic Martin

Catherine Droitcourt

PhD Laurence Gusdorf

Claire Abasq

Ali Dadban

Christine Chiaverini

Pierre Vabres

PhD Denis Herbreteau

MD Olivia Boccara

Michel Wassef

MD Annabel Maruani
email: annabel.maruani@univ-tours.fr

Annabel Pr

Maruani

Intramuscular capillary-type hemangioma: Diagnosis, treatment, and outcomes. A French multicentric retrospective study of 66 cases

Keywords: Intramuscular capillary-type hemangioma: diagnosis, treatment, and outcomes Intramuscular capillary-type hemangioma: diagnosis, treatment, and outcomes Intramuscular hemangioma, Hemangioma of the skeletal muscle, Intramuscular capillary-type hemangioma, Intramuscular fast-flow vascular anomaly, Intramuscular vascular malformation, Extracranial arteriovenous malformation

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

INTRODUCTION

Intramuscular capillary-type hemangiomas (ICTHs) belong to the group of intramuscular "hemangiomas", a misnomer often used to name different types of intramuscular vascular malformations. ICTH has been reported under several names including intramuscular hemangioma (or angioma) small-vessel type or capillary-type and intramuscular infiltrating angiolipoma [START_REF] Allen | Hemangioma of skeletal muscle, An analysis of 89 cases[END_REF][START_REF] Pribyl | 2 nd , Infiltrating angiolipoma or intramuscular hemangioma? A report of five cases[END_REF][START_REF] Beham | Intramuscular angioma: a clinicopathological analysis of 74 cases[END_REF][START_REF] Yilmaz | Intramuscular capillary-type hemangioma: radiologic-pathologic correlation[END_REF]. Its first description in 1972 by Allen and Enzinger, based on pathology, divided "hemangiomas of skeletal muscle" into 3 subgroups according to the predominant vessel size in the lesion: the small-vessel, large-vessel, and mixed types [START_REF] Allen | Hemangioma of skeletal muscle, An analysis of 89 cases[END_REF]. The small-vessel type, or ICTH, was then better identified with a recent description of specific MRI features [START_REF] Yilmaz | Intramuscular capillary-type hemangioma: radiologic-pathologic correlation[END_REF][START_REF] Orly | Characteristics, natural course and treatment of intramuscular capillarytype hemangioma: A systematic literature review[END_REF], but its diagnosis remains challenging.

ICTH is a benign vascular anomaly of the skeletal muscle, consisting of an intramuscular solid mass that most often grows gradually, without pain, and mostly affects young people. Its evolution is poorly understood, but it does not regress spontaneously as does infantile hemangioma [START_REF] Yilmaz | Intramuscular capillary-type hemangioma: radiologic-pathologic correlation[END_REF][START_REF] Orly | Characteristics, natural course and treatment of intramuscular capillarytype hemangioma: A systematic literature review[END_REF]. The term "hemangioma" for ICTH might be confusing, because it is used for infantile hemangioma and congenital hemangioma, 2 very different conditions of newborns and infants [START_REF] Wassef | ISSVA Board and Scientific Committee, Vascular anomalies classification: Recommendations from the International Society for the Study of Vascular Anomalies[END_REF]. The main differential diagnoses of ICTH are other intramuscular vascular lesions, mainly intramuscular venous malformations, phosphatase and tensin homolog (PTEN)-related hamartomas of soft tissue (PHOSTs), and fibro-adipose vascular anomalies (FAVAs). These correspond to the large vessels and most of the mixed subtype hemangiomas of the skeletal muscles as described by Allen and Enzinger [START_REF] Allen | Hemangioma of skeletal muscle, An analysis of 89 cases[END_REF]. Venous malformations, often called "cavernous hemangiomas" or "venous hemangiomas" are slowflow malformations that can be located in the muscles, often containing phleboliths, which induce pain. On MRI, they present as well-lobulated lesions with low signal intensity (SI) on T1-weighted images (WIs) and high SI on T2-WIs and delayed post-contrast images enhancement on T1-gadolinium sequences. Communication with the normal deep venous system throughout narrow tributaries can be seen. In some cases, thrombosed venous pouches or phleboliths may mimic signal voids [START_REF] Hein | Venous malformations of skeletal muscle[END_REF][START_REF] Eryilmaz | Intramuscular cavernous hemangioma of the temporalis muscle[END_REF]. PHOSTs are recently described vascular entities that often present as diffuse infiltrative intramuscular tumors referred to as PTEN hamartomas that may contain slow-flow and/or fast-flow vascular abnormal vessels. They might be isolated or associated with various systemic features of PTEN hamartoma tumor syndrome [START_REF] Kurek | PTEN hamartoma of soft tissue: a distinctive lesion in PTEN syndromes[END_REF]. Radiologically, PHOSTS show a less clearly defined mass with high SI on T2-WIs and large flow voids due to arterial supply. The lesion might potentially involve surrounding tissues. More recently, FAVAs were described; these are additional rare vascular entities characterized by an infiltrating solid lesion that contains slow-flow vessels and a fatty component that usually involves muscular and neurovascular structures and is associated with severe pain and joint contracture [START_REF] Alomari | Fibro-adipose vascular anomaly: clinical-radiologic-pathologic features of a newly delineated disorder of the extremity[END_REF][START_REF] Lipede | Management of Fibro-adipose Vascular Anomalies (FAVA) in paediatric practice[END_REF]. On MRI, they can be recognized on T1-WIs, showing features of fatty and fibrous tissue components with striated multilayer features similar to fatty muscle infiltration. On T2-WIs, a hyper SI with "patchy spots" might be observed, with serpiginous structures (venous and lymphatic component) with low SI T2-WI foci (thrombi or phleboliths). Enhancement is moderate to marked after gadolinium injection.

ICTH is a fast-flow vascular entity (which is a key feature for diagnosis), with currently no gold standard for diagnosis. Somatic mutations of the MAP2K1 and KRAS genes, also found in extracranial arteriovenous malformations (AVMs), were identified in a series of ICTH cases, which raises the question of the distinction between ITCH and intramuscular AVMs [START_REF] Goss | Intramuscular fast-flow vascular anomaly contains somatic MAP2K1 and KRAS mutations[END_REF]. However, these molecular results still need to be confirmed in larger cohorts. The International Society for the Study of Vascular Anomalies, which classified ICTHs as "provisionally unclassified vascular anomalies" [START_REF] Wassef | ISSVA Board and Scientific Committee, Vascular anomalies classification: Recommendations from the International Society for the Study of Vascular Anomalies[END_REF] without specifying whether they belong to vascular tumors or vascular malformations, considered the possible overlap between ICTH and intramuscular AVMs, both belonging to a spectrum of fast-flow intramuscular lesions.

In this study, we aimed to collect all cases of ICTH followed up in 9 tertiary care centers for vascular anomalies in France and to assess the diagnostic criteria for better classification, treatments and outcomes.

MATERIAL AND METHODS

Design and setting

This was a multicenter retrospective observational study of all ICTH cases followed up from January 1, 2000 to March 1, 2021 in France. The study was approved by the ethics committee of Tours (#2021 014).

Participants and samples

We included all patients with a lesion diagnosed by their hospital radiologist as an intramuscular hemangioma, with at least MRI images (T1, T2, T1-injected sequences) available. Because of no current gold standard for diagnosing ICTH, an adjudication group of experts was constituted to establish the diagnosis by taking into account a range of arguments (imaging and pathological data). All eligible cases were reviewed by the adjudication committee, consisting of a pathologist (MW), 2 radiologists (AB, GB), 1 surgeon (AJ) and 2 dermatologists (OB, AM) who reviewed US reports, MRI images, CT-scan and angiogram images when available and pathological reports. Criteria for inclusion in the final cohort of ICTH patients were that all experts agreed with the diagnosis of ICTH on the basis of data compatible with their knowledge and experience of the entity and agreed with one another. In case of disagreement or no consensus, the case was excluded, as were cases for which another diagnosis was established.

Outcomes

History, clinical, radiological, pathological and molecular data were analyzed.

Diagnostic criteria

Clinical characteristics that fit the diagnosis of ICTH were the occurrence of an intramuscular tissular mass without extension to the skin. Cases presenting skin necrosis, a port-wine stain, a bluish skin discoloration, a segmental overgrowth or any kind of hamartoma were excluded.

On MRI, lesions should be intramuscular lesions delineated or not by a fatty border, hyperintense on T2-weighted images and isointense to skeletal muscle on T1-weighted images, with or without an intralesional fatty component. On T1-weighted sequences, one should observe hypervascularization and homogeneous contrast enhancement after gadolinium injection. Cases have no calcifications (phleboliths), which suggests a venous malformation. When ultrasonography is performed, criteria are a non-compressible mass within the skeletal muscle, iso-or echogenic, and color and spectral Doppler ultrasonography showing arterial-type vascularization. When performed, angiography could show enlarged arterial feeder arteries with opacification of the tissue mass before any contrast appears in the draining veins and no arteriovenous nidus. Pathological features that fit with a ICTH diagnosis were: presence of a proliferation of capillary-type vessels separating or infiltrating between the skeletal muscles, and absence of clusters of large ectatic vessels, calcifications (phleboliths), and bundles of smooth-muscle cells. Molecular data were not mandatory, although they were collected.

Statistical analyses

Statistical analyses were mainly descriptive. Some data were missing because of the retrospective nature of the study. Categorical variables are expressed with number (percentages) and quantitative variables with median (interquartile range [IQR]). When comparison tests were performed, we considered p<0.05 as statistically significant. We used R v3.6.2 for analysis.

RESULTS

Among 133 patients assessed for eligibility, with history obtained from 9 French centers, 66 were included (Figure 1, flow chart). The main reasons for exclusion were insufficient data collection to allow for a consensual diagnosis (n=23), non-intramuscular lesion (n=10), and other diagnoses (n=34).

Table 1 summarizes patient characteristics. The median age at the first signs was 28.0 years (IQR 21.036.0; range: birth to 58 years). The male/female ratio was 1.28. The median time to diagnosis was 12 months (IQR 5.022.0).

ICTH was most commonly a progressive increasing mass (83.9%), very rarely warm (3.1%), and with a firm consistency (83.3%); it was mostly painless (88.9%) (Table 1). A triggering factor was suggested in 8 (12.1%) cases, in particular pregnancy (n=4) and local trauma (n=3). ICTH could occur all over the body, but the most affected body regions were the head and neck area (n=28/66, 42.4%). Overall, among the 61 cases that included a description of the affected muscle(s), 31 muscles or muscle groups were affected (Figure 2). The most commonly affected muscle was the masseter (n=12).

Imaging features are reported in Table 2. On MRI, the most common features were a well-delineated mass (87.8%) that was homogeneous (60.5%), isointense on T1-weighted images (63.6%), always displaying hypersignals on T2-weighted images, and after contrast enhancement, showing high and homogenous contrast enhancement (92.3%). Most cases (83.3%) showed flow voids. Angiography was performed in 13 (19.7%) cases, mostly describing a tissular blush with no nidus, inhomogeneous parenchymal staining without early draining or arteriovenous shunting. Ultrasonography was performed in 46 (69.7%) cases, describing mainly a hyperechogenic intramuscular lesion (87.5%) with a high resistance index (55.6%). Biopsy or surgical excision was performed in 51 (77.3%) cases, was not performed in 12 cases and was unknown for 3 (Table 1). Descriptions mainly involved capillary proliferation (46.7%), with small-size and concomitant other-sized vessels in 86.1% of cases and presence of adipose tissue in 75.7% (Figure 3). A lobulated architecture was reported in 35.5% of cases and endothelial cells between muscle fibres in 17.9%. On immunostaining, when described, all cases were negative for GLUT-1 and positive for ERG, AML, CD31 and CD34. The Ki67 proliferation index was < 10% in all cases. In one case, somatic mutations of the MAP2K1 and KRAS genes was searched, with negative results.

Among the 66 ICTH cases, the consensus group identified 59 cases of typical ICTH and 7 atypical cases that showed features overlapping AVM and ICTH, which we called AVM-like ICTH. The latter were more voluminous than typical ICTH cases and could involve several muscles. They were located on the neck (n=1), trunk (n=2) and limbs (n=4) but not on the head, contrary to typical ICTH for which the face location was frequent (Table 3). The pathological description of both types was similar. MRI features of AVM-like ICTH were close to those of typical ICTH except that the lesions were larger, less well-delineated, and more heterogeneous, with constant flow voids in our cases and moderate to high contrast enhancement (Figures 4 &5). Angiography revealed high vascularization with early venous return as in AVM and no nidus but global parenchymal vascularization such as in typical ICTH, with a feeding dysplasic artery and mild arteriovenous shunting.

The median post-treatment follow-up was 17.0 months (IQR 8.5-31.0). The main therapeutical procedure was surgical resection: 11/47 (23.4%) patients with available data underwent complete surgery alone and 6/47 (12.8%) underwent complete surgery after embolization, which was curative in all cases, with no recurrence (Table 4). Embolization alone (n=6, 12.8%) resulted in complete regression in one case. In 13/47 (27.7%) cases, clinical/imaging follow-up showed no spontaneous regression. Sirolimus, prescribed at 0.1 mg/day for 6 months in one patient was inefficient. Three serious adverse events occurred after treatment: major bleeding that required transfusion during partial surgical removal alone, bilateral cerebellar ischemia after embolization alone of a lesion located in the neck, and haemorrhagic shock during surgery the day after embolization. No patient died from disease, explorations or treatment.

DISCUSSION

This retrospective multicenter French cohort of ICTH is the largest reported in the literature. The 66 cases included were mostly young adults and the lesions presented as a progressively enlarging firm and painless mass, mostly affecting the head and neck region (especially the masseter muscle) and limbs. Seven cases shared some imaging features with AVM (less welldelimited tissue mass, larger tortuous afferent arteries, more early draining-vein opacification, mild arteriovenous shunting) and could represent an intermediate lesion between ICTH and AVM. Clinically, these latter lesions tended to be larger, less frequently located on the face, and more painful and barely accessible for complete surgery. We propose to name them AVM-like ICTH. As compared with AVM-like ICTH, extracranial AVM might be apparent at birth; the lesions increase in size progressively, sometimes more suddenly during hormonal changes, and tend to recur after surgery [START_REF] Liu | extracranial arteriovenous malformations: natural progression and recurrence after treatment[END_REF]. The presence of rapid arteriovenous shunting, with early opacification of the veins before any tissue opacification, enlarged tortuous arterial feeders and multiple collateral supply are the main discriminating features for a positive diagnosis of extracranial AVM [START_REF] Geibprasert | Imaging characteristics of arteriovenous malformations simulating vascular tumors of the temporalis muscle: a report of two cases[END_REF]. Goss et al., identifying same KRAS and MAPK somatic mutations in ICTH as in extracranial AVMs, need to be confirmed [START_REF] Goss | Intramuscular fast-flow vascular anomaly contains somatic MAP2K1 and KRAS mutations[END_REF], because they would have potential impact, beyond the nosological aspect, on the future use of targeted therapies for non-operable mutated ICTH [START_REF] Nicholson | Successful management of an arteriovenous malformation with trametinib in a patient with capillary-malformation arteriovenous malformation syndrome and cardiac compromise[END_REF]. Only one case was investigated for these mutations in our cohort, with negative results. Further molecular investigations are needed of both typical and AVM-like ICTHs.

In the literature, intramuscular hemangioma is misnamed for actually describing different intramuscular lesions. In this cohort, about half of the cases screened were excluded because the committee of experts reclassified them as other entities, the most frequent being intramuscular venous malformations. These lesions are vascular anomalies that appear as heterogeneous masses on MRI, are not fast-flowing, and can contain phleboliths, a special type of organized thromboses [START_REF] Hein | Venous malformations of skeletal muscle[END_REF][START_REF] Fowell | Venous malformations of the head and neck: current concepts in management[END_REF][START_REF] Seront | Venous malformations of the head and neck[END_REF][START_REF] Dompmartin | Association of localized intravascular coagulopathy with venous malformations[END_REF]. Usually, they appear as well-delineated spherical hypointense areas on T1-and T2-weighted MRI images, which argues for a venous malformation. FAVAs and PHOSTs were other differential diagnoses. In FAVA, MRI reveals a heterogeneous lesion containing rounded areas with high signal spots associated with fibroadipose components, corresponding to large dysplastic veins, typical of this disease. PHOST often presents as irregularly shaped intramuscular lesions, showing T1-weighted contrast enhancement and T2-hyperintense signals, with dilated, serpiginous vascular components and small arteriovenous fistulae with dilated draining veins [START_REF] Kurek | PTEN hamartoma of soft tissue: a distinctive lesion in PTEN syndromes[END_REF][START_REF] Alomari | Fibro-adipose vascular anomaly: clinical-radiologic-pathologic features of a newly delineated disorder of the extremity[END_REF][START_REF] Sheybani | Fat-containing soft-tissue masses in children[END_REF][START_REF] Arnold | Diagnostic imaging of vascular anomalies[END_REF][START_REF] Pattamapaspong | Imaging of intramuscular haemangiomas of the extremities[END_REF]. Both FAVA and PHOST correspond histologically to lesions described as "angiomatosis of softtissue" by Rao and Weiss [START_REF] Rao | Angiomatosis of soft tissue. An analysis of the histologic features and clinical outcome in 51 cases[END_REF].

ICTHs are benign lesions that mainly grow fast until becoming stable in growth. Eight of our 66 cases featured a potential triggering factor such as a trauma and/or pregnancy. Four additional cases reported in the literature occurred after pregnancy [START_REF] Orly | Characteristics, natural course and treatment of intramuscular capillarytype hemangioma: A systematic literature review[END_REF]. This observation might be fortuitous but needs to be kept in mind because trauma and hormonal changes are also potential triggering factors for extracranial AVMs [START_REF] Liu | extracranial arteriovenous malformations: natural progression and recurrence after treatment[END_REF][START_REF] Azis | Extracranial arteriovenous malformations rupture in pregnancy[END_REF]. Complete surgical resection, when possible, preceded or not by embolization, seems curative, with very low recurrence risk.

This study has some limitations, such as biases linked to its retrospective nature, which inevitably induces missing data.

In conclusion, clinical history and physical examination, associated with MRI images, allow for the diagnosis of ICTH in most cases. The most important characteristics of these lesions are that they are delineated masses presenting intralesional flow voids and early homogenous T1 gadolinium enhancement. In rare cases, a biopsy is required to exclude malignant soft tissue tumors. We identified a subgroup of large ICTH lesions displaying more heterogeneous MRI features that we called AVM-like ICTH. These lesions require angiography to distinguish them from extracranial AVMs.

TABLES Table 1. Clinical, pathological and molecular data for intramuscular capillary hemangioma (ICTH) (n=66 patients)

  signs, years, median [IQR] Age at diagnosis, years, median [IQR] Diagnostic delay, months, median [IQR]

Figure 1 .Figure 2 .Figure 3 .

 123 FIGURE LEGENDS

Figure 4 .

 4 Figure 4. Typical intramuscular capillary hemangioma of the scalene muscle in a 25-year-old patient: (A) Firm intramuscular mass underlying a normal skin;(B-C) MRI showing a welllimited homogeneous T2-weighted hyperintense signal (B) on a coronal T2 sequence with fat suppression and strong and homogenous enhancement following gadolinium injection (C) as seen on a T1-weighted-sequence with fat suppression. On a digital subtraction angiography with bolus injection in the ascending cervical artery, the lesion is seen as an early blush (D) demonstrating its richly vascularized nature; at the intermediate phase (E), the lesion shows a pseudo-mesh appearance with vascular thickening along the main axis of muscular fibers.

Figure 5 .

 5 Figure 5. Arteriovenous malformation (AVM)-like intramuscular capillary hemangioma: Clinically evident as a (A) large firm mass of the scapulae area underlying a normal skin; (B-D) MRI images showing the lesion with a T1 signal isointense to that of adjacent paravertebral muscle (B); T2 hyperintense lobulated lesion within the paravertebral muscles (C, arrowheads) and heterogenous intramuscular infiltrating enlargement with multiple hypointense pseudo-tubular areas, corresponding to high flow vessels (flow voids, C, black circles). (D): A frontal projection of a digital subtraction angiography following left subclavian artery injection shows an early arterial intralesional blush; in an ascending branch of the costo-cervical trunk (D, inset) a (possibly flow-related) arterial aneurysm (white arrow) upstream of the lesion parenchyma.

Figure 1 -

 1 Figure 1. Flow chart

Table 2 . Imaging data for ICTH (n=66 patients)

 2

	Variable	Description
	MRI, n (%)	66 (100)
	Description of the mass
	No consensus on delineation	17
	Well-delineated mass	43 (87.8)
	Heterogeneous mass	17 (39.5)
	Homogeneous mass	26 (60.5)
	No consensus on homogeneity	23
	T1-weighted sequence description (n=22)
	Missing data	44
	Hyperintense signals	3 (13.6)
	Hypointense signals	5 (22.7)
	Isointense signal	14 (63.6)
	T2-weighted sequence description (n=58)
	Missing data	8
	Hyperintense signals	58 (100)
	Contrast enhancement (n=39)
	Missing data	27
	High contrast	36 (92.3)
	High and homogeneous contrast	13
	Moderate contrast	3 (7.7)
	Flow voids (n=54)	45 (83.3)
	Missing data	12
	CT scan, n (%)	15 (22.7)
	Angiography, n (%)	13 (19.7)

Table 3 . Comparison of clinical and imaging data between typical intramuscular capillary-type hemangioma (ICTH) and AVM-like ICTH

 3

		Typical ICTH	AVM-like ICTH	Missing	p
		(n=59)	(n=7)	Data
	PATIENT CHARACTERISTICS
	Sex, n (%)			0	-
	Females	29 (43.9)	3 (42.9)
	Males	37 (56.1)	4 (57.1)
	Age at first signs, years, median	28.17 [22.3, 36]	19.0 [6.75, 33.5]	13	0.361
	[IQR]
		30 [23.9, 38.5]	21 [11.1, 31.2]	0	0.148
	Age at diagnosis, years, median
	[IQR]	12.0 [6.0, 22.5]	6.0 [1.2, 16.0]	13	0.231
	Diagnostic delay, months, median
	[IQR]
	CLINICAL DATA
	Lesion size			37	0.464
	Area (mm 2 ), median [IQR]	1950.0 [1080.5-	3675.0 [2812.5-
		4212.5]	4537.5]
	Pain, n (%)	4 (7.1)	3 (42.9)	3	0.025
	Reported discomfort, n (%)	7 (12.3)	2 (33.3)	3	0.201
	Topography, n (%)			0	0.350
	Head and neck	27 (45.8)	1 (14.3)
	Trunk	8 (13.6)	2 (28.6)
	Lower limb	12 (20.3)	2 (28.6)
	Upper limb	12 (20.3)	2 (28.6)
	Progressive increase, n (%)	48 (87.3)	4 (57.1)	29	0.062

Table 4 . Management and outcomes of ICTH (n=47 patients)

 4

		ICTH,
	Outcome, n (%)	intramu
		scular
		capillar
		y
		hemang
		ioma

DECLARATION OF INTERESTS

This research received no external financial or non-financial support.

There are no additional relationships to disclose.

There are no patents do disclose.

There are non additional activities to disclose.


domain: sdv math

text:
Remy Ben Messaoud

Vincent Le Du

Brigitte Charlotte Kaufmann

Baptiste Couvy-Duchesne

Lara Migliaccio

Paolo Bartolomeo

Mario Chavez

Fabrizio De

Vico Fallani

Brigitte Charlotte

Baptiste Couvy- Duchesne

Low-dimensional controllability of brain networks

Keywords: Network controllability, Spectral graph theory, Brain networks, Diffusion weighted imaging, Hemispheric lateralization

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

Introduction

Network controllability refers to the possibility of identifying the driver nodes in an interconnected system and opportunely modifying their dynamics to steer the system's state towards desired configurations 1 . Network controllability stems from other established disciplines, such as control and graph theory, and has recently witnessed significant theoretic advances. Analytical approaches allowed to establish the impact of network topology on controllability by unveiling the link between the minimum number of drivers necessary to control a network and its node degree distribution 2,3 as well as its eigenvalues' geometric multiplicity [START_REF] Yuan | Exact controllability of complex networks[END_REF] . More recently, the use of controllability metrics derived from Gramian-based energy minimization provided analytical tools to measure the potential of candidate drivers to physically control networks [START_REF] Pasqualetti | Controllability Metrics, Limitations and Algorithms for Complex Networks[END_REF][START_REF] Wang | Physical controllability of complex networks[END_REF][START_REF] Lindmark | Minimum energy control for complex networks[END_REF] . The number of applications studying the controllability of real systems has flourished in the last years. By modeling how information is stipulated to pass along the underlying network structure, controllability informs on the dynamical properties of the nodes that cannot be merely obtained by looking at their connectivity alone. Structural controllability can quantify the ability to manipulate and regulate the flow in power grids [START_REF] Yang | Critical Nodes Identification of Complex Power Systems Based on Electric Cactus Structure[END_REF] or transportation networks [START_REF] Rinaldi | Controllability of transportation networks[END_REF] , influence information and behavior in social networks [START_REF] Cremonini | Controllability of social networks and the strategic use of random information[END_REF] , and gain insights into complex biological processes, such as gene expression regulation [START_REF] Rajapakse | Dynamics and control of state-dependent networks for probing genomic organization[END_REF] and brain functions [START_REF] Yan | Network control principles predict neuron function in the Caenorhabditis elegans connectome[END_REF][START_REF] Gu | Controllability of structural brain networks[END_REF][START_REF] Tang | Colloquium: Control of dynamics in brain networks[END_REF] . Despite its great potential, network controllability's concrete impact remains quite limited due to the difficulty of identifying in a reliably not only the number of drivers but more importantly where they are located and what is the control signal that is to be injected to push the system into desired states. This is a typical computationally expensive problem, further aggravated by the extremely large condition numbers that arise when the driver nodes are less than one-fifth of the entire network [START_REF] Sun | Controllability Transition and Nonlocality in Network Control[END_REF] . In those very common situations, round-off numerical errors arise and the network cannot be controlled in practice despite placing the drivers strategically to ensure theoretical controllability [START_REF] Wang | Physical controllability of complex networks[END_REF][START_REF] Sun | Controllability Transition and Nonlocality in Network Control[END_REF][START_REF] Yan | Spectrum of controlling and observing complex networks[END_REF] . Notably in neuroscience, studying the control influences between differently specialized brain systems is not only important from a fundamental perspective (endogenous control) [START_REF] Bressler | Large-scale brain networks in cognition: emerging methods and principles[END_REF] but often a necessity imposed by the existing technology that allows to stimulate only one or few regions at a time (exogenous control) [START_REF] Dayan | Noninvasive brain stimulation: from physiology to network dynamics and back[END_REF][START_REF] Momi | Network-level macroscale structural connectivity predicts propagation of transcranial magnetic stimulation[END_REF] . However, because of the aforementioned problems, the typical controllability metrics are often close to zero or negative, making it difficult to predict the actual impact on brain functioning and resulting behavior [START_REF] Tu | Warnings and Caveats in Brain Controllability[END_REF][START_REF] Suweis | Brain controllability: Not a slam dunk yet[END_REF][START_REF] Stocker | Network controllability measures of subnetworks: implications for neurosciences[END_REF] . To address these limitations disparate approaches have been introduced. From an algorithmic perspective, using a stepwise exploration to measure the extent of the network that can be controlled by a driver, provides some benefit on the overall accuracy [START_REF] Bassignana | Stepwise target controllability identifies dysregulations of macrophage networks in multiple sclerosis[END_REF] . Targeting specific parts of the system is perhaps the most intuitive way to lower the imbalance between the number of drivers and the size of the network to control [START_REF] Gao | Target control of complex networks[END_REF][START_REF] Klickstein | Energy scaling of targeted optimal control of complex networks[END_REF] . In general, output controllability can be adopted to reduce the problem dimensionality by controlling condensed components of the original network state. For example, controlling the averaged state of the nodes within separate groups, or clusters, has been demonstrated to be an effective strategy to improve precision [START_REF] Casadei | Model Reduction Based Approximation of the Output Controllability Gramian in Large-Scale Networks[END_REF] . However, average-state controllability is an over-approximation of the real dynamics as nodes in the same cluster might have very different functional states. Because network dynamics are constrained by the underlying structure [START_REF] Boccaletti | Complex networks: Structure and dynamics[END_REF][START_REF] Porter | Dynamical Systems on Networks: A Tutorial[END_REF] , one possibility would be instead to exploit the existing connectivity to have a more representative aggregation. By leveraging spectral graph theory, we introduced a framework that embeds the original network state into few representative components characterizing the structure of the system at different organizational scales [START_REF] Chung | Spectral graph theory[END_REF] . We formally introduced our mathematical approach for linear-time invariant systems and evaluated via extensive simulations on synthetic networks the control precision and representativeness in contrast to the standard approaches. Eventually, by projecting the Gramian into a low-dimensional space we introduced a parsimonious controllability metric to identify the driver nodes of a network. We used such low-dimensional controllability metric to identify the most influential areas of the human connectome obtained by gathering more than 6000 samples from the UK-biobank cohort [START_REF] Sudlow | UK Biobank: An Open Access Resource for Identifying the Causes of a Wide Range of Complex Diseases of Middle and Old Age[END_REF] and provide a map of the causal interactions between functional brain systems and characterize cerebral specialization.

Results

Without loss of generality, we considered unweighted and undirected networks whose dynamics are governed by the following linear time-invariant equation

"(%) = ("(%) + *+(%) ,(%) = -"(%) (1)
Where " ∈ ℝ ! contains the states of the 0 nodes, ( ∈ ℝ !×! is a stabilized version of the network adjacency matrix 1 and * ∈ ℝ !×! ! directs the external input + ∈ ℝ ! ! into 0 $ ≤ 0 distinct driver nodes. Here, , ∈ ℝ % with 3 < 0 represents the low-dimensional output state obtained by projecting the entire state vector " into fewer components via the so-called output matrix -∈ ℝ %×!

,(%) = > % " =(%) = -)*+ "(%) (2)
where the output matrix -)*+ = > % 8 ( is obtained by selecting and reordering a number 3 < 0 of eigenmaps via a filtering matrix > % (Fig 1b).

Tradeoff between control precision and representativeness

We considered hierarchical modular small-world (HMSW) networks whose Laplacian eigenmaps informed on the network partition across multiple scales [START_REF] Sporns | Small-world connectivity, motif composition, and complexity of fractal neuronal connections[END_REF]

(Fig 2a, Materials and Methods).

Without loss of generality [START_REF] Luenberger | Introduction to dynamic systems: theory, models, and applications[END_REF] , we set the initial network state in the origin " , = 0 @⃗ ∈ ℝ ! and the final state as a random departure from it " -~C(D -= 1 , F - ' = 100) ∈ ℝ ! , which corresponded to steering the related eigenstate from , , = 0 @⃗ to , -= -)*+ " -. Candidate drivers were chosen by ranking the nodes according to their betweenness centrality so as to ensure homogeneous distances from the rest of the network [START_REF] Duan | Prevalence and scalable control of localized networks[END_REF] . We then computed the associated input signals +(%) by minimizing the output cost function

G . (+, % -) =∥ , --,(% -) ∥ ' + I ∫ ∥ +(K) ∥ ' LK / " ,
, where I is a regularization parameter weighting the accuracy of the solution with respect to the signal energy [START_REF] Shirin | Optimal control of complex networks: Balancing accuracy and energy of the control action[END_REF] (Text S1). By injecting the optimal inputs back into the model dynamics, we considered the cosine similarity [START_REF] Senoussaoui | A Study of the Cosine Distance-Based Mean Shift for Telephone Speech Diarization[END_REF] to measure the control accuracy in terms of precision M = 0(/ " ) • 0 " ∥0(/ " )∥∥0 " ∥ . Here, I = 10 56 and % -= 1 ensured sufficiently accurate solutions (Fig S1a). To evaluate the impact of the low-dimensional control on the original network state, we also measured control accuracy in terms of representativeness N = 7(/ " ) • 7 " ∥7(/ " )∥∥7 " ∥ .

Results showed that the number of eigenmaps significantly affects the control precision. The lower was 3, the higher was M regardless of the number of drivers 0 $ . However, having more drivers led to higher precision and faster transitions towards the largest values (Fig 2b). To simulate more realistic situations, we considered the case of single-driver control. On average, reducing the number of eigenmaps improved the control precision but worsened representativeness. Optimal tradeoffs could be then explored considering unbiased indicators such as the total accuracy M + N (Fig 2c). We next targeted the control of a subnetwork O with P ≤ 0 nodes and computed the Laplacian taking into account its internal connectivity (Fig 2d, Materials and Methods). As in full network control, we computed the related low-dimensional eigenstates and measured the total accuracy as a function of the target size. Results showed an overall benefit in controlling smaller parts of the network, especially when considering internal drivers and using a relatively low number of eigenmaps (Fig 2e).

Note that these findings were obtained by selectively including eigenmaps from the most to the least informative in terms of , magnitude. Changing the inclusion criterion, such as ranking them according to the eigenvalues, led to similar tendencies albeit with less representativeness (Fig S1b). Notably, results stayed comparable even when considering networks with different topologies, i.e., random and scale-free (Fig S1c).

Low-dimensional controllability of brain systems

To quantify the ability of individual brain regions to influence the activity of other areas, we introduced a control metric based on a low-dimensional projection of the Gramian matrix 5

Q = R S 89 ** ( S 8 # 9 LK : , (3)
where * selects the i th candidate driver. Similar to Eq. 2, we operated such projection via the -)*+ matrix, i.e., Q )*+ = -)*+ Q-)*+ ( . Its smallest eigenvalue T ;<! )*+ gave a low-dimensional version of the worst-case control centrality T ;<! , measuring the amount of energy needed for node U to steer the network into the most difficult-to-reach state [START_REF] Baggio | Energy-Aware Controllability of Complex Networks[END_REF] .

We considered structurally weighted brain networks obtained from healthy human subjects in the UK-biobank cohort. Brain networks consisted of 0=214 nodes including cortical and subcortical regions of the Schaefer atlas 39 (File S1). They were grouped into nine differently specialized functional systems [START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF] (Fig 3a, Materials and Methods). Our first main result showed that reducing the number of eigenmaps allowed the entire brain to be controllable from every single node (i.e., T ;<! )*+ > 0) whereas this failed when using standard T ;<! . The transition to full controllability occurred at 3 * ≤ 5 in line with previous evidence [START_REF] Sun | Controllability Transition and Nonlocality in Network Control[END_REF] and used in all our subsequent analyses (Fig 3b). Compared to standard worst-case controllability some regions exhibited similar relative spatial contributions (e.g. SUB), while others faded out or gained importance notably in the visual-temporal part of the right hemisphere (e.g.

VIS, DAN) (Fig 3c).

Using the target controllability approach described in the previous section, we next considered the Laplacian of the subnetworks corresponding to the brain systems and evaluated their controllability (Fig S2). Results pooled from all subjects showed that systems could be controlled from both internal and external nodes enabling to derive optimal driver-target configurations (Fig 4a, Tab 1). More in general, we reported a system controllability 〈T ;<! )*+ 〉 hierarchy with subcortical (SUB) and visual (VIS) ones among the easiest to control (Fig 4b). Internal drivers exhibited higher control centrality compared to external ones, indicating a preferential system selfregulation 〈T ;<! )*+ 〉 <! with respect to external regulation 〈T ;<! )*+ 〉 >?/ (Materials and methods). That was particularly evident for SUB, but also for the primary systems VIS and somatomotor (SMN), the latter to a lesser extent (Fig 4c). This result was not necessarily due to the spatial proximity between the driver and the target, but rather to the topological distance in terms of shortest paths (Fig S2).

We then investigated the system lateralization in terms of self-regulation to accumulate in one or the other hemisphere. For each targeted system, we computed a lateralization index Z = @ $ 5@ % @ $ A@ % , where [ = 〈T ;<! )*+ 〉 <! for the nodes in the right (R) or left (L) hemisphere. Results showed significant system-dependent lateralization, which was not merely due to the possible differences in terms of system size between hemispheres (Fig 4d). Such lateralization was perfectly in line with the known functional experimental evidence, such as the right-dominated limbic and dorsal-attention network, and could not be retrieved when using the standard worstcase control metric T ;<! .

Mapping inter-system controllability in the brain

To better understand how different systems were influencing each other we next studied their reciprocal role as drivers and targets. For each pair of systems i and j we established a directed weighted link by calculating the geometric mean control centrality of the nodes in i when targeting j, i.e., 〈T ;<! )*+ 〉 <→C . By inspecting the resulting meta-graph, we found a heterogeneous strength and distribution of causal influences that enabled to identify the preferential targets and drivers for each system (Fig 5a). Notably, primary systems exerted a strong controlling influence over the dorsal-attentional network (VIS®DAN, SMN®DAN), while more balanced interactions were found between associative systems (e.g., FPCN«DMN, DMN«LIM).

To quantify the control unbalance, we computed for each system the difference between the outgoing and incoming sum of the weighted links of the meta-graph. Results confirmed that both VIS and SMN exhibited a significantly higher tendency to control as compared to SUB and DAN which instead appeared much easier to be controlled (Fig. 4b). A more balanced profile was instead observed for the other systems, such as the limbic system (LIM). Similar results were also obtained when using a finer-grained version of the Yeo2011 parcellation (Fig. S4). We finally investigated hemispheric dominance in terms of control capacity. To do so, we computed the T ;<! )*+ value of each node when targeting the two hemispheres, separately. We defined i) the ipsilateral control as the mean 〈T ;<! )*+ 〉 <→< of the nodes in one hemisphere targeting the same hemisphere, and ii) the contralateral control as the mean 〈T ;<! )*+ 〉 <→C of the same nodes targeting the other hemisphere.

Taken separately, each hemisphere exhibited a stronger capacity to control itself in contrast to its ability to influence the other. However, non-trivial results emerged when comparing the two hemispheres. While the ipsilateral control of the left hemisphere (LàL) was significantly higher than the right hemisphere (RàR), the contralateral control of the right hemisphere (RàL) was significantly higher than the left hemisphere (LàR) (Fig 5c). This global behavior could be explained by the same tendency across all the local brain systems, but the temporoparietal junction (TPJ) whose contralateral control exhibited an opposite trend (Fig 5d). Altogether, these findings indicated that the left hemisphere exhibits better self-regulation, while the right hemisphere has a stronger ability to influence its counterpart.

Discussion

Spectral graph theory and network controllability Spectral graph theory (SGT) is a powerful and versatile framework for understanding and analyzing complex interconnected systems. SGT is particularly concerned with the eigenvalues and eigenvectors of the Laplacian matrix associated with a graph and informs on the global and finer-grained connectivity structures in the network [START_REF] Villegas | Laplacian renormalization group for heterogeneous networks[END_REF][START_REF] Gfeller | Spectral Coarse Graining and Synchronization in Oscillator Networks[END_REF] . Laplacian eigenvectors, also known as eigenmaps, find various applications in network science, such as data clustering and dimensionality reduction via spectral embedding [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] . From a broader graph signal processing angle, the Laplacian eigenmaps technique informs on dynamical processes like synchronization or diffusion and is a suitable framework to analyze how external perturbations propagate in a network [START_REF] Boccaletti | Complex networks: Structure and dynamics[END_REF] . More recently, SGT has been developed to identify connectivity gradients in the brain [START_REF] Huntenburg | Large-Scale Gradients in Human Cortical Organization[END_REF] , disentangle brain structure and function [START_REF] Preti | Decoupling of brain function from structure reveals regional behavioral specialization in humans[END_REF] , improve machine learning for systems neuroscience [START_REF] Richiardi | Machine Learning with Brain Graphs: Predictive Modeling Approaches for Functional Imaging in Systems Neuroscience[END_REF] , and decompose brain dynamics into graph Fourier modes [START_REF] Lioi | Gradients of connectivity as graph Fourier bases of brain activity[END_REF][START_REF] Huang | A Graph Signal Processing Perspective on Functional Brain Imaging[END_REF] .

Despite its potential to reduce the complexity of signals on graphs, SGT has been poorly explored in the context of network controllability. Network controllability provides a mathematically grounded framework to identify drivers, influence network dynamics and understand causal relationships, but it still faces critical challenges that hamper its usability in many practical situations [START_REF] Pasqualetti | Controllability Metrics, Limitations and Algorithms for Complex Networks[END_REF][START_REF] Sun | Controllability Transition and Nonlocality in Network Control[END_REF][START_REF] Yan | Spectrum of controlling and observing complex networks[END_REF] . This is mainly due to the presence of ill-posed conditions that make the problem numerically unstable with important consequences on the solutions, such as an overestimation of the drivers or completely unreliable input-controlling signals [START_REF] Sun | Controllability Transition and Nonlocality in Network Control[END_REF] . Hence, the development of methods to reduce the problem complexity is paramount to enhance the usability and interpretation of network controllability in real-world applications. Targeting specific subnetworks [START_REF] Bassignana | Stepwise target controllability identifies dysregulations of macrophage networks in multiple sclerosis[END_REF][START_REF] Gao | Target control of complex networks[END_REF][START_REF] Klickstein | Energy scaling of targeted optimal control of complex networks[END_REF] or reducing the network to fewer averaged nodes [START_REF] Casadei | Model Reduction Based Approximation of the Output Controllability Gramian in Large-Scale Networks[END_REF][START_REF] Souza | Controlling complex networks with complex nodes[END_REF] have shown substantial benefits, yet they represent approximations that oversimplify the original network structure.

SGT offers a more flexible solution to reduce problem complexity enabling to write the network states as a linear combination of its actual connectivity structures (the eigenmaps) across different topological scales. By establishing a link between SGT and output controllability, we showed that controlling the most representative spectral components from the graph Fourier transform, instead of the original network state, leads to significant improvements even in the presence of unbalanced \ 0 $ P ] ≪ 1_ and noisy configurations. The price to pay for such a control precision consisted of an almost equivalent loss in the representativeness of the solution in terms of capacity to reproduce the original network state. Optimal trade-offs can be further investigated in future research to better capture the interplay between network topology, the number of output components, and the smoothness of the state solution.

Theoretical and effective controllability of brain networks

Network controllability can provide insights into how different regions of the brain interact to perform various cognitive and motor functions by passing electrical activity through white matter connections [START_REF] Tang | Colloquium: Control of dynamics in brain networks[END_REF][START_REF] Honey | Predicting human resting-state functional connectivity from structural connectivity[END_REF] . Its theoretical framework can be used to experimentally test basic questions in network neuroscience and inform noninvasive brain stimulation or neurofeedback approaches to enhance cognitive functions, such as memory and learning [START_REF] Cornblath | Temporal sequences of brain activity at rest are constrained by white matter structure and modulated by cognitive demands[END_REF] . Furthermore, investigating the controllability of brain networks can shed light on the causal mechanisms underlying neurological [START_REF] Zhou | Predicting regional neurodegeneration from the healthy brain functional connectome[END_REF] and psychiatric disorders [START_REF] Braun | Brain network dynamics during working memory are modulated by dopamine and diminished in schizophrenia[END_REF] as well as help identify potential entry points for therapeutic interventions, such as brain stimulation techniques for conditions like epilepsy, Parkinson's disease, stroke and chronic pain 1 .

Despite its potential, brain network controllability faces a number of methodological challenges that limit its usability. This is mainly due to the unreliability of control centrality metrics, such as the worst-case control T ;<! , which are mathematically correct but have little translational value when the number of drivers is much lower than the nodes to control. The main consequence is that is not always possible to infer the control power of single brain areas because it cannot be assessed numerically and remains difficult to interpret [START_REF] Tu | Warnings and Caveats in Brain Controllability[END_REF][START_REF] Suweis | Brain controllability: Not a slam dunk yet[END_REF][START_REF] Stocker | Network controllability measures of subnetworks: implications for neurosciences[END_REF] . Here, we showed that lowdimensional controllability allows overcoming this limitation and reveals previously unappreciated brain drivers in the VIS and DAN that could not be identified by merely looking at their connectivity strength (Fig 3,Fig S5). SUB was the most controllable system possibly due to its rich-club organization integrating information from remote parts of the brain [START_REF] Bell | Subcortical contributions to large-scale network communication[END_REF] , and together with VIS and SMN exhibited a high tendency to self-regulate as compared to other associative systems [START_REF] Menardi | Maximizing brain networks engagement via individualized connectome-wide target search[END_REF] .

The analysis of the self-regulation's lateralization allowed to specify the known functional [START_REF] Corbetta | Control of goal-directed and stimulus-driven attention in the brain[END_REF] and anatomical [START_REF] De Schotten | A lateralized brain network for visuospatial attention[END_REF] hemispheric asymmetries such as the right preference for TPJ and SVAN partially overlapping with the ventral attention network, and for LIM in line with the control of emotional processing [START_REF] Gainotti | Emotions and the Right Side of the Brain[END_REF] . In addition, the found left-hemisphere advantage of SMN and DMN aligned with the suggested preferential role in action selection [START_REF] Schluter | Cerebral dominance for action in the human brain: the selection of actions[END_REF] and endogenously generated cognition [START_REF] Liu | Hemispheric asymmetries in visual mental imagery[END_REF] . More importantly, we provided new insights for those systems whose asymmetry is more controversial such as VIS [START_REF] Liu | Evidence from intrinsic activity that asymmetry of the human brain is controlled by multiple factors[END_REF] and DAN [START_REF] Shulman | Right Hemisphere Dominance during Spatial Selective Attention and Target Detection Occurs Outside the Dorsal Frontoparietal Network[END_REF][START_REF] Nobre | The attentive homunculus: Now you see it, now you don't[END_REF] . From a control-theoretic perspective, these systems exhibited a strong tendency to self-regulate via the right hemisphere and provided a novel anatomical basis for the common clinical observation that spatial neglect is more frequent, persistent, and severe after damage to the right hemisphere as compared to lesions to the left homotopic areas [START_REF] Lunven | Attention and spatial cognition: Neural and anatomical substrates of visual neglect[END_REF] .

Causal relationships between brain systems

The presence of several systems, often called "networks", coexisting in the brain is an emerging feature of its intrinsic functioning and crucial to understanding human behavior in both healthy and pathological conditions [START_REF] Corbetta | Control of goal-directed and stimulus-driven attention in the brain[END_REF] . By leveraging connectomics to trace functional connectivity, several large-scale brain systems have been identified that are important, among others, for decision making (FPCN), emotion, motivation and memory (LIM), attention for externallydirected tasks (DAN), introspective thoughts (DMN), visual and somatosensory processing (VIS, SMN), and internal/external thinking (SVAN) [START_REF] Baker | A Connectomic Atlas of the Human Cerebrum-Chapter 1: Introduction, Methods, and Significance[END_REF] . Notably, these systems exhibit a hierarchy to their operation, but they do not work in isolation. Instead, they are known to integrate and synchronize to carry out complex functions [START_REF] Bressler | Large-scale brain networks in cognition: emerging methods and principles[END_REF][START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF] .

At a larger scale, the brain hemispheres can also be considered as anatomically separated systems that functionally interact to accomplish complex tasks [START_REF] Gazzaniga | Cerebral specialization and interhemispheric communication: Does the corpus callosum enable the human condition?[END_REF][START_REF] Doron | Dynamic network structure of interhemispheric coordination[END_REF][START_REF] Martínez | Role of interhemispheric connections in functional brain networks[END_REF] . Cerebral specialization is not only important to underpin the neural basis of language, attention, and motor control, but it is also critical to assess and treat related deficits due to mental or neurological disorders [START_REF] Gazzaniga | Cerebral specialization and interhemispheric communication: Does the corpus callosum enable the human condition?[END_REF][START_REF] Bartolomeo | From competition to cooperation: Visual neglect across the hemispheres[END_REF] . Yet, determining how large-scale systems simultaneously influence each other in a causal manner is still poorly understood. Here, we showed that drivers in the primary systems (VIS, SMN) were particularly apt to influence distributed activity in associative and attentional systems (e.g., DAN, SVAN), while drivers in the prefrontal cortex exhibited facilitated control of SUB, SMN, and DMN (Fig 4). This result was further confirmed by the meta-wiring diagram showing how brain systems simultaneously interact from a control theoretic perspective.

The role of primary systems as controlling drivers and that of associative systems as controlled targets, match previous results obtained with network communication models of brain cognition [START_REF] Srivastava | Models of communication and control for brain networks: distinctions, convergence, and future outlook[END_REF] and align with theories postulating that the brain topography is organized along spatial and functional gradients [START_REF] Margulies | Situating the default-mode network along a principal gradient of macroscale cortical organization[END_REF][START_REF] Malkinson | From perception to action: Intracortical recordings reveal cortical gradients of human exogenous attention[END_REF] . These gradients span from peripheral regions that handle perception and action to core areas dealing with more abstract cognitive functions. Our findings revealed that core attentional networks are in particular controlled by primary peripheral systems, which are instead more influenced by external events. In terms of hemispheric asymmetries, our results extend recent evidence using fMRI functional connectivity suggesting a preference for the left hemisphere to interact with itself and for the right hemisphere to integrate information both locally and from the other hemisphere [START_REF] Gotts | Two distinct forms of functional lateralization in the human brain[END_REF] . Our findings provide an anatomical basis for this fundamental result by unveiling unbalanced directed influences within and between hemispheres.

Conclusion

Network controllability is not only crucial to identify driver nodes, but in general to explore causal relationships in complex interconnected systems such as the human brain. By leveraging principles from spectral graph theory, the introduced low-dimensional controllability framework offers a reliable tool to interrogate, influence and ultimately understand the behavior of complex systems.

Materials & Methods

Laplacian of stabilized networks. The graph Laplacian 5 = 6 -( is a symmetric matrix that can be rewritten via singular value decomposition as 5 = 8 ( Λ8, where Λ contains its increasing eigenvalues (T & ≤ T ' … ≤ T ! ) and 8 = [8 & , 8 ' , … , 8 ! ] are the associated eigenvectors8 [START_REF] Chung | Spectral graph theory[END_REF] . By construction, eigenvectors are orthonormal (i.e., 8 ( 8 = a) and constitute a projection basis to obtain the spectral representation of any signal " on a graph via the so-called Graph Fourier Transform (GFT) [START_REF] Chung | Spectral graph theory[END_REF] , i.e., " = = 8 ( ". Here, networks corresponded to graphs whose adjacency matrices 1 coded for the presence/weight of links between the nodes. To ensure the stability of the system and the existence of the Gramian we next added negative self-loops to the network via a linear transformation ( = 1 -ba, with the constant b chosen here as the maximum eigenvalue of 1 ( " E (%) and derived its low-dimensional output

, E )*+ (%) = > % " E e (%) = -E )*+ " E (%) (4)
where the output matrix -E )*+ = > % 8 E ( is obtained by selecting and reordering a number 3 < P of spectral components via a filtering matrix > % .

Hierarchical modular small-world network and controllability metrics. We used the HMSW model [START_REF] Sporns | Small-world connectivity, motif composition, and complexity of fractal neuronal connections[END_REF] to generate synthetic networks with 0 = 256 nodes, 8 = log ' (256) hierarchical levels, and an average network density of 0.035. Other parameters were: initial cluster size=2, and connection density fall-off per level=2.5. The choice of this model and its parameters was particularly relevant for the purposes of this study and has been widely adopted in topological analysis of brain networks [START_REF] Sporns | Small-world connectivity, motif composition, and complexity of fractal neuronal connections[END_REF] which exhibit hierarchical modularity [START_REF] Meunier | Modular and Hierarchically Modular Organization of Brain Networks[END_REF] .

Starting from the worst-case low-dimensional control centrality T ;<! )*+ of each node, we defined aggregated metrics to measure larger topological effects, i.e., i) system controllability 〈T ;<! )*+ 〉 as the mean of all the nodes targeting a specific subnetwork, ii) self-regulation 〈T ;<! )*+ 〉 <! as the mean of all the nodes inside a target subnetwork, and iii) external regulation 〈T ;<! )*+ 〉 >?/ as the mean of all the nodes outside a target subnetwork.

Neuroimaging data and tractography. The UK-biobank dataset is a large-scale cohort containing clinical, genetic, and imaging data [START_REF] Sudlow | UK Biobank: An Open Access Resource for Identifying the Causes of a Wide Range of Complex Diseases of Middle and Old Age[END_REF] . We selected N=6134 human subjects who had both T1 weighted and diffusion MRI and no known disease history, i.e., international classification of diseases ICD-10='none'. Population statistics: 50.73% women, 88.72% righthanded, and age 62.41±7.25 at the time of the MRI scans. Informed consent was obtained from all UK Biobank participants. We rejected those who requested the withdrawal of their data. Procedures are controlled by a dedicated Ethics and Guidance Council (http://www.ukbiobank.ac.uk/ethics), with the Ethics and Governance Framework available at https://www.ukbiobank.ac.uk/media/0xsbmfmw/egf.pdf. IRB approval was also obtained from the North West Multi-centre Research Ethics Committee. This research has been conducted using the UK Biobank Resource under Application Number 53185. Downloaded imaging data were corrected and preprocessed with the UK-biobank pipeline [START_REF] Alfaro-Almagro | Image processing and Quality Control for the first 10,000 brain imaging datasets from UK Biobank[END_REF] .

Tractography and connectome construction. As for the processing, we computed Tissue response and Fiber Orientation Distribution using multi-tissue and multi-shell algorithms in MRtrix3 80 . T1 images were aligned to an extracted mean b0 volume via the FLIRT function in FSL 81 . We performed a 5-tissue type segmentation to compute the grey-white matter interface. These were then used in MRtrix3 to compute an anatomically constrained tractography with a cut-off of 0.1 and a density of 1M streamline that was shown to be sufficient for reproducibility [START_REF] Roine | Reproducibility and intercorrelation of graph theoretical measures in structural brain connectivity networks[END_REF] . As for the parcellation, we used the Schaefer atlas of 200 cortical regions (i.e., nodes) exhibiting high structural and functional representativeness [START_REF] Luppi | Combining network topology and information theory to construct representative brain networks[END_REF] . Its regions were grouped into eight cortical systems according to the Yeo2011 mapping 40 : VIS, SMN, DAN, SVAN, LIM, FPCN, DMN, TPJ. The cortical atlas was complemented by 14 subcortical regions from the FreeSurfer segmentation that form the network SUB [START_REF] Menardi | Maximizing brain networks engagement via individualized connectome-wide target search[END_REF] . The parcellation was transferred to the subject space using the T1 linear co-registration and the UK-biobank warp field. It was finally dilated and masked to be used in MRtrix3 along the SIFT [START_REF] Smith | The effects of SIFT on the reproducibility and biological accuracy of the structural connectome[END_REF] for the connectome extraction. Fiber assignment was done with a radial search of 3mm and the resulting connectomes were symmetric with a zero diagonal. Link weights correspond to the number of fibers between two nodes.

Reproducibility

All the simulations for both synthetic and brain networks were done in MATLAB and the code is made available at github.com/Inria-NERV/Network-Control. The HMSW model and shortest paths were computed with the BCT toolbox (bctnet.com) [START_REF] Rubinov | Complex network measures of brain connectivity: Uses and interpretations[END_REF] . All neuroimaging data are available upon request from the UK-biobank (ukbiobank.ac.uk) and the code of the processing pipeline is available at github.com/UKB-dwi-2-connectome.

Figure 1 -Principles of low-dimensional network controllability

a) Toy network with 0=7 nodes. Its state " = [" &, " ' , … , " F ] can be seen as a signal over a graph. By exploiting the Laplacian eigenvectors 8 = [8 & , 8 ' , … , 8 F ] of the network, the signal " can be embedded in a spectral space via the graph Fourier transform (GFT) " l = 8 ( ". The resulting eigenstate " l measures how the signal is spatially distributed across different topological scales, from coarser (" l & , " l ' ,… ) to finergrained (…, " l G , " l F ) ones. We can also define the inverse operation iGFT " = 8" l that is illustrated in this panel. b) Low-dimensional controllability in terms of linear-time invariant (LTI) network output control. Instead of focusing on the network state, the goal is to determine the input signal +(t) that steers a low-dimensional output given by the subset of the eigenstate ,(%) = > % " =(%) = -)*+ "(%), where the output matrix -)*+ = > % 8 ( is obtained by selecting and reordering a number 3 < 0 of spectral components via a filtering matrix > % . In this example, the network has 0=7 nodes and the first 3 = 3 spectral components are selected and reordered arbitrarily. Here, the filtering matrix > % is a 3 × 7 matrix whose elements are pℎ &,' r = pℎ ',& r = pℎ H,H r = 1 and zero elsewhere.

b)

a)

= ! " ! × ! % ! % " % # 1 -1 0 + ! " " × + ! " # × … Driver Output ' ( ) $%& = + ' % ( ,(.) 0 ! . = ! " " (.) 0 " (.) = ! " ! (.) 0 ) (.) = ! " ) (.)
• Stationary equation:

0 = 𝜕ℋ 𝜕𝑢 = 2𝜌𝑢(𝑡) + 𝐵 𝑇 𝑣(𝑡) (4)
• Boundary/Transversality condition: 𝑣(𝑡 𝑓 ) = 𝜕𝐸(𝑥(𝑡 𝑓 ))

𝜕𝑥(𝑡 𝑓 )

= -2𝐶 𝑇 ( 𝑦 𝑓 -𝐶𝑥(𝑡 𝑓 ))

The adjoint and stationary equations, (3) and (4), can be rewritten in condensed form:

[ 𝑥 * v * ̇] = [ 𝐴 -(2𝜌) -1 𝐵𝐵 𝑇 0 -𝐴 𝑇 ] [ 𝑥 * 𝑣 * ] = 𝐻 [ 𝑥 * 𝑣 * ] ( 6
)
where 𝐻 is the Hamiltonian matrix. Equation ( 6) can be solved as:

This way, we obtained an expression for 𝑣 * (𝑡) and the problem is almost solved since equation

The first two equations come from (7) and the third is the Hamiltonian boundary condition (5).

We have three unknowns 𝑣 * (0), 𝑥 * (𝑡 𝑓 ), 𝑣 * (𝑡 𝑓 ), and three equations. The calculus gives:

𝑣 * (0) = (𝑀 22 -2𝐶 𝑇 𝐶𝑀 12 ) † ((2𝐶 𝑇 𝐶𝑀 11 -𝑀 21 )𝑥 0 -2𝐶 𝑇 𝑦 𝑓 ), (9) where † denotes the pseudo-inverse of a matrix.

Laplacian of the network and state matrix

For a graph with adjacency matrix 𝐺 the Laplacian matrix is defined as 𝐿(𝐺) = 𝐷 𝐺 -𝐺, where 𝐷 𝐺 is the diagonal matrix containing the degree sequence of nodes of 𝐺. In this work, we stabilize the nodes' dynamics and consider the matrix 𝐴 = 𝐺 -𝑐𝐼, where 𝑐 is an arbitrary constant. We verify here, using the linearity of the Laplacian operator, that 𝐿(𝐴) = 𝐿(𝐺). b) Adjacency matrix of the group-averaged controllability meta-graph between the different brain systems. The links represented by color are the group-average of control centrality of the system in column i when targeting system in line j 〈𝜆 𝑚𝑖𝑛 𝐸𝐼𝐺 〉 𝑖→𝑗 . 〈𝜆 𝑚𝑖𝑛 𝐸𝐼𝐺 〉 𝑖→𝑗 is obtained by taking the geometric mean control centrality 𝜆 𝑚𝑖𝑛 𝐸𝐼𝐺 of drivers in system i when targeting system in j. Self-loops and the SUB network are not represented as their control centrality is several orders of magnitudes higher.

Figure 2 -Figure 3 -Figure 4 -Figure 5 -

 2345 Figure 2 -Fundamental advantage of low-dimensional controllability a) Laplacian eigenvectors (or eigenmaps) for one realization of the hierarchical modular small-world network (HSWM). Colors identify the different spatial contributions of the nodes at increasingly finer topological scales. b) Control accuracy in terms of precision (M) as a function of the number of eigenmaps. Different lines correspond to different number of drivers. Candidate drivers were chosen by ranking the nodes according to their betweenness centrality. Results correspond to the mean values obtained from 100 simulated networks. c) Average node precision ⟨M⟩ = & ' ∑ M < ! <I&

( 4 )

 4 gives 𝑢 * (𝑡) = -𝐵 𝑇 𝑣(𝑡) 2𝜌 . We still are left with the unknown 𝑣 * (0), which can be found using the boundary conditions3 . If we note 𝑒 𝑡 𝑓 𝐻 = [ 𝑀 11 𝑀 12 𝑀 21 𝑀 22 ], we have the three following boundary constraints: { 𝑥 * (𝑡 𝑓 ) = 𝑀 11 𝑥 0 + 𝑀 12 𝑣 * (0) 𝑣 * (𝑡 𝑓 ) = 𝑀 21 𝑥 0 + 𝑀 22 𝑣 * (0) 𝑣 * (𝑡 𝑓 ) = -2𝐶 𝑇 ( 𝑦 𝑓 -𝐶𝑥 * (𝑡 𝑓 ))

  Fig S4 -Control relationships between brain systems in a finer parcellation. a) The 17 split components of the Yeo2011 brain atlas parcellation.

  c) System total contribution as the sum of outgoing and incoming weighted links from the individual meta-graphs. Bars indicate group-average values and error bars standard deviations. d) System control unbalance as the difference between the sum of outgoing and incoming weighted links from the individual meta-graphs. Positive values=tendency to act as driver. Negative value= tendency to act as target. Bars indicate group-average values and error bars standard deviations.

  Fig S5 -Relationship between low-dimensional controllability and node strength. a) Group-averaged spatial distribution of node strength k given by the sum of all the weighted links of a node. b) Group-averaged spatial distribution of low-dimensional (𝜆 𝑚𝑖𝑛 𝐸𝐼𝐺 ) control centrality. c) Scatter plot of node strength k and low-dimensional (𝜆 𝑚𝑖𝑛 𝐸𝐼𝐺 ) control centrality for all 𝑛 × 𝑁 nodes and participants. Pearson correlation test revealed a low correlation between the two metrics.

  [START_REF] Anderson | Optimal Control: Linear Quadratic Methods[END_REF][START_REF] Wang | Controllability of networked MIMO systems[END_REF] . In general, different configurations for the output matrix can be chosen to aggregate the original network state into fewer components. Instead of selectingin an arbitrary or predetermined fashion, we derived it directly from the network structure by considering the Laplacian linear operator 5 = 6 -(, where 6 is a diagonal matrix containing the node degree sequence[START_REF] Chung | Spectral graph theory[END_REF] . By construction, 5 can be decomposed into several eigenvectors 8 = [8 & , 8 ' , … , 8 ! ] or eigenmaps, associated with incrementing eigenvalues that inform on the network structure load from coarser to finer-grained topological scales. Eigenmaps serve as a basis to derive a spectral representation of the network state " = = 8

(

", here referred to as eigenstate (Fig

1a)

. By selecting a smaller number of eigenmaps, we aimed to control the low-dimensional eigenstate

  E is a diagonal matrix that contains the degree sequence of the target nodes excluding the links to the rest of the network. Note that a precise relation holds with the submatrix of the original Laplacian 5(O, O) = 5 E + d E , where d E is a diagonal matrix containing the sum of the links connecting the target nodes to the rest of the network 77 . Hence, 5 E exactly matches the actual spatial topological scales of the subnetwork and it is not biased by external influences. By considering the state of the subnetwork " E (%) and the related eigenmaps 8 E we obtained the target eigenstate " E

	e (%) =
	8 E

[START_REF] Anderson | Optimal Control: Linear Quadratic Methods[END_REF][START_REF] Karrer | A practical guide to methodological considerations in the controllability of structural brain networks[END_REF]

. It can be easily demonstrated that such transformation does not alter the network Laplacian which is a linear operator (Text S1).

Target network control. Let

O = [c & , c ',

… , c

; ] be a subset of P ≤ 0 target nodes indexed by c < ∈ [1, 0] and 1(O, O) the corresponding adjacency matrix. We formulated the related target control problem by considering the Laplacian of the subnetwork 5 E = 6 E -1(O, O), where 6

Acknowledgments

This work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement N° 864729). BCD is supported by a CJ Martin Fellowship (NHMRC app 1161356). The authors declare no competing interests.

Target system

Inside driver Outside driver

Table 1: Driver-target best configurations.

Drivers are selected according to their maximum group-averaged low-dimensional control centrality T ;<! )*+ .

Supplementary text S1

Optimal control: input signal derivation In the results section, we validated our framework for synthetic networks by simulating trajectories. We first, need to find the control input 𝑢(𝑡). To do so we solve the following optimization problem: The problem can be solved using Pontriyargin's maximum principle 1 (see sections 5.1 and 5.2) by introducing the Hamiltonian equation:

where 𝑣(𝑡) is the vector of the adjoint states. We know from the fundamentals of optimal control theory 2 that the optimal trajectory (𝑥 * , 𝑣 * , 𝑢 * ) is the solution to the following equations: b) Representativeness 𝜂 and precision 𝛿 (inset) as a function of the final state dispersion 𝜎 𝑓 using 𝑟=64 most representative eigenmaps sorted according to the magnitude of the eigenstate (in blue) or using the 𝑟 first eigenmaps corresponding to the smaller Laplacian eigenvalues (in red). We used all nodes as drivers 𝑛 𝑑 =𝑛 and the final state had a fixed mean of 𝜇 𝑓 =1. for smooth final states with small small dispersion, both eigenmaps selection schemes gave high and similar representativeness 𝜂 with a slight superiority of the 𝜆-sorted scheme. As the dispersion increased, 𝜎 𝑓 >1, selecting the eigenmaps accordingly to the eigenstate gave better representativeness 𝜂 and that was accompanied with lower precision 𝛿 in the low-diensional space as the control task was harder. Results are averaged over 100 HMSW network realizations and shaded areas represent standard deviations. c) Precision 𝛿, representativeness 𝜂, and their sum as a function of the number of controlled eigenmaps 𝑟 for two different topologies: the Erdos-Renyi model ER(𝑛=256 , 𝑝=0.035), and the Barabasi-Albert model BA(𝑛=256 and connection density 0.035, bias 𝛾= 2). d de. We also the same control task as for the HMSW model: the final state was sampled from 𝒩(𝜇 𝑓 =1, 𝜎 𝑓 2 =100), and the trajectory was simulated for 𝜌=10 -4 and 𝑡 𝑓 =1.

Supplementary figures (S1-S5)

Results are averaged over 100 network realizations.


Hung Duong

Vassili Pham

Denis Pustovalov

Kouamé

The Performance Improvement of Ultrasound Localization Microscopy (ULM) Using the Robust Principal Component Analysis (RPCA)

Keywords: ULM, ultrafast ultrasound, inverse problems, robust PCA, blood flow, tissue filtering, UCA localization

This paper presents an algorithm for ultrafast ultrasound localization microscopy (ULM) used for the detection, localization, accumulation, and rendering of intravenously injected ultrasound contrast agents (UCAs) enabling to yield hemodynamic maps of the brain microvasculature. It consists in integrating a robust principal component analysis (RPCA)based approach into the ULM process for more robust tissue filtering, resulting in more accurate ULM images. Numerical experiments conducted on an in vivo rat brain perfusion dataset demonstrate the efficiency of the proposed approach compared to the most widely used state-of-the-art method.

I. INTRODUCTION

Ultrasound (US) imaging is increasingly used in medical practice as a non-invasive, nonionizing radiation, low-cost, and ease-of-use tool. Using plane wave techniques, the US modality can perform acquisitions with a frame rate up to 20kHz, which enables a good visualization and quantification of microvascular structure including sudden functional changes in the brain, as well as an increase of blood-low sensitivity for more accurate filtering [START_REF] Tanter | Ultrafast imaging in biomedical ultrasound[END_REF]. However, the propagation of US waves faces many phenomena including diffraction and attenuation, resulting in a limited spatial resolution. The recent advent of ultrasound localization microscopy (ULM) successfully handles the trade-off of US imaging between spatial resolution and penetration depth while achieving a high sensitivity thanks to the clever interaction between sparse ultrasound contrast agents (UCAs) and ultrafast imaging [START_REF] Couture | Ultrasound Localization Microscopy and Super-Resolution: A State of the Art[END_REF]. Having a similar principle to photoactivated localization microscopy (PALM), ULM considers UCAs as sensors to image the region of interest by detecting, localizing, accumulating, and rendering thousands of microbubble events through density-based techniques. Besides, the high framerate within ULM allows the retrieval of the trajectories of UCAs in the bloodstream, thus imaging efficiently both the vasculature and microvasculature. Interestingly, ULM has been freshly adapted for 2D transcranial application, which has the potential of providing new imaging biomarkers for early detection, diagnosis, and prognosis of many cerebral diseases [START_REF] Demené | Transcranial ultrafast ultrasound localization microscopy of brain vasculature in patients[END_REF]. Despite its efficiency, most current algorithms for ULM have been implemented together with a trivial spatio-temporal singular value decomposition (SVD) filter [START_REF] Demené | Transcranial ultrafast ultrasound localization microscopy of brain vasculature in patients[END_REF]- [START_REF] Heiles | Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy[END_REF]. The idea of using this technique emanates from the fact that in US imaging, particularly ultrafast technique, D. H. Pham, V. Pustovalov and D. Kouamé are with the IRIT Lab, Université de Toulouse, and CNRS UMR 5505, 31400 Toulouse, France. Email: {duong-hung.pham; vassili.pustovalov; denis.kouame}@irit.fr.

back-scattered blood signals carrying UCAs (originating from vessels) and clutter signals (from tissues surrounding the vessels) superimpose themselves on the images [START_REF] Demené | Spatiotemporal Clutter Filtering of Ultrafast Ultrasound Data Highly Increases Doppler and fUltrasound Sensitivity[END_REF]. Applying tissue filtering consequently enables the removal of unwanted clutter signals from stationary or moving tissues to retain only the contribution of individual UCAs. The SVD technique is fast and straightforward; however, the empirical tuning process of its thresholding parameters does not rely on any interpretable model, which remains problematic in practice, limiting the method's applicability.

In this paper, we propose to use an alternative technique for tissue filtering, namely a robust principal component analysis (RPCA), resulting in more accurate and saturated ULM images [START_REF] Candès | Robust principal component analysis?[END_REF]. The principle of RPCA is to formulate tissue filtering as an inverse problem of the blood-tissue separation, where a priori knowledge about tissue and blood structures carrying UCAs can be incorporated. In practice, the tissues are assumed to be low-rank while the blood flow is sparse. The inverse problem is then solved using iterative approaches such as the alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. The remainder of the paper is organized as follows. The background about SVD-based ULM is given in Section II. The proposed PRCA-based ULM algorithm is detailed in Section III. Finally, numerical results on in vivo ultrafast US data are regrouped in Section IV, showing the improvement achieved by the proposed approach over the most efficient algorithm for ULM nowadays.

II. BACKGROUND ON SVD-BASED ULM

Ultrasound localization microscopy (ULM) is a multi-step process including six main steps resumed as follows 1 . 1) Acquisition of raw US images of a perfused organ injected with UCAs at a high frame rate. 2) Tissue filtering on these images to get rid of tissue signals. 3) Detection and localization of UCAs with sub-pixel precision. 4) Pairing successive positions of the detected UCAs together into tracks. 5) Interpolation and projection of the tracks on a rendering grid. 6) Reconstruction of ULM vascular image accumulating a large number of tracks (called rendering). The quality of ULM images is strongly dependent upon all the steps above which have been greatly improved over the past decade. For instance, a thorough study of ULM framework has been introduced including 6 sets of in vitro and in vivo data and a systematic comparison of the seven most used algorithms for UCAs localization in [START_REF] Heiles | Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy[END_REF]. In this article, we focus on the improvement of the tissue filtering step which is often overlooked in previous works [START_REF] Demené | Transcranial ultrafast ultrasound localization microscopy of brain vasculature in patients[END_REF]- [START_REF] Heiles | Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy[END_REF]. Indeed, in these works, a trivial SVD was used to separate blood flow containing UCAs from unwanted tissue signals, and also to increase the Doppler sensitivity of the resulting image. It consists in factorizing the 2D Casorati matrix S ∈ C NzNx×Nt (i.e. obtained by stacking the frames into columns of the recorded image) using SVD as follows [START_REF] Demené | Spatiotemporal Clutter Filtering of Ultrafast Ultrasound Data Highly Increases Doppler and fUltrasound Sensitivity[END_REF]:

S = U ΣV † = r k=1 u k σ k v † k , (1)
where

N z , N x ,
B = T b k=Tc u k σ k v † k . (2)
Although the SVD technique is very simple and can perform at a very low computational cost, the manual choice of the optimal thresholds is not obvious, not only in the ULM scenario but also in other practical situations, thus making its applicability complex.

III. PROPOSED RPCA-BASED ULM

The main aim of the proposed algorithm is to enhance the ultimate ULM image quality by integrating a more efficient tissue filtering technique in the ULM process in the place of SVD. To this end, we investigate a robust principal component analysis (RPCA)-based technique [START_REF] Candès | Robust principal component analysis?[END_REF] as an alternative solution. It consists in formulating blood-tissue separation as an inverse problem, where a priori knowledge about tissue and blood structures with UCAs is taken into account. The blood flow with UCAs B is assumed to be sparse and is conventionally promoted by the minimization of the l 1 -norm. The tissue signals T is assumed to change very slowly with time, leading to a low-ranking assumption on T which is usually modeled by the nuclear norm ||.|| * . N is assumed to be a Gaussian noise. Thus, RPCA is expressed as the following convex optimization problem:

[ B, T ] = arg min B,T ||S -B -T || 2 F + λ||B|| 1 + ρ||T || * , (3)
where . stands for the estimated variables, ||.|| F is the Frobenius norm, and λ, ρ > 0 are two hyperparameters balancing the trade-off between the blood sparsity and the tissues' low-rankness [START_REF] Wright | Robust Principal Component Analysis: Exact Recovery of Corrupted Low-Rank Matrices via Convex Optimization[END_REF]. To solve (3), the augmented Lagrangianbased alternating direction method of multipliers (ADMM) is used. It consists in solving iteratively several sub-problems over each variable separately [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. Precisely, the augmented Lagrangian related to (3) is first written as follows:

L (B, T , ν) = λ||B|| 1 + ρ||T || * + µ 2 ||S -B -T + 1 µ ν|| 2 F ,
where ν is the Lagrange multiplier and µ is the Lagrangian penalty parameter controlling the convergence speed of the algorithm. Finally, at each iteration k, ADMM performs the following three steps, until a predefined stopping criteria is met:

B(k+1) = arg min B (λ||B||1 + µ 2 ||B -(S -T (k) + 1 µ ν (k) )|| 2 F ) T (k+1) = arg min T (ρ||T || * + µ 2 ||T -(S -B (k+1) + 1 µ ν (k) )|| 2 F ) ν (k+1) = ν k + µ(S -B (k+1) -T (k+1) ).
It is interesting to note that the two first steps for B(k+1) and T (k+1) are both convex problems associated, respectively, with closed-form solutions: soft thresholding (ST) [START_REF] Lin | The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices[END_REF] and singular value thresholding (SVT) [START_REF] Cai | A Singular Value Thresholding Algorithm for Matrix Completion[END_REF]. Note that other RPCA-related approaches, e.g. [START_REF] Pham | Joint Blind Deconvolution and Robust Principal Component Analysis for Blood Flow Estimation in Medical Ultrasound Imaging[END_REF], [START_REF] Pham | Fast High Resolution Blood Flow Estimation and Clutter Rejection Via an Alternating Optimization Problem[END_REF] can also be used.

IV. NUMERICAL RESULTS

This section regroups numerical results on a in vivo rat brain perfusion dataset to illustrate the improvement brought by the proposed algorithm over the most efficient method currently for ULM. Note that due to space limitation, only this in-vivo data will be used to illustrate the proposed approach. As mentioned above, an exhaustive study conducted in [START_REF] Heiles | Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy[END_REF] drew a conclusion that SVD-based UML with radial symmetry (RS) algorithm for UCAs localization came out on top both in qualitative and quantitative evaluations. Therefore, it will be selected for comparison with our proposed algorithm in the sequel.

A. In vivo rat brain perfusion dataset

This dataset is the in vivo angiography of the rat brain that is publicly available at [START_REF] Arthur | OPULM PALA, Type: dataset[END_REF]. Before the acquisition, the skull of the rat had been removed following a craniotomy surgery protocol. Then, UCAs (400µl, Bracco Imaging, Italy) were continuously injected through the jugular vein (192, 000 US images). The radiofrequency data were beamformed on the ultrasound machine's graphics processing unit (SuperSonic Imagine). An US acquisition comprising of 240 blocs of 800 frames of size 78x118 pixels was taken at a compounded frame rate of 1KHz with 3 tilted plane waves [-3 o , 0 o , +3 o ] and pulse repetition frequency 15 MHz (Vermon, France) with a 0.1mm pitch. The acquisition time was 4 min. Refer to [START_REF] Heiles | Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy[END_REF] for more details.

B. Performance comparison

The hyperparameters of the two studied algorithms were tuned to maximize the performance measures. To this end, those of SVD-based ULM were kept the same as in the original publication, which showed the best performance [START_REF] Heiles | Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy[END_REF].

  and N t are respectively the axial, lateral, and temporal dimensions. U ∈ C NzNx×NzNx and V ∈ C Nt×Nt are two unitary matrices of the spatial (left) u k and temporal (right) v k singular vectors. Σ ∈ R NzNx×Nt is a non-square diagonal matrix with diagonal elements being its singular values σ k . † denotes the conjugate transpose and r = min(N z N x , N t ) the rank of S. It is worth mentioning that the diagonal elements Σ are sorted in ascending order of magnitude. The smallest singular values correspond to the noises, the intermediate ones to the blood signals, and the largest ones to the tissue signals. Then, some empirical heuristics are used to tune the clutter T c and blood T

b thresholds. This enables to filter out the unwanted tissue signals and the noises, and also to reveal of the blood flow components as follows:

The interested reader is referred to[START_REF] Heiles | Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy[END_REF] for a detailed description of these steps of the method.

Particularly, the SVD relying on two distinct thresholds to filter the eigenvalues were, respectively, chosen as T c = 5 and T b = 800. For the RPCA-based ULM, we recall that λ and ρ associated with RPCA express the compromise between the blood sparsity and the tissue low-rankness while µ manages the convergence speed of the algorithm. As suggested in [START_REF] Wright | Robust Principal Component Analysis: Exact Recovery of Corrupted Low-Rank Matrices via Convex Optimization[END_REF] for the general RPCA problem, ρ was set to 1 while the reference values

and µ ref = 10 × λ ref , were applied to ease the tuning. With this in mind, the best values of λ and µ selected were respectively set to 0.0104 and 0.104, in the experiments. Note also that the RS algorithm for UCAs localization with the same parameter setting for both two studied algorithms was used to ensure a fair comparison.

In Fig. 1, we depict the ULM density rendering images along with two zoomed-in regions obtained by the two algorithms carried out on the in vivo rat brain perfusion dataset. The density rendering image shows the number of trajectories passing through a single pixel whose square size is set to 10µm×10µm. Visually judging these plots, one may notice that the proposed RPCA-based ULM algorithm yields a density rendering image that is overall less noisy and more saturated vessels than the SVD-based ULM one. Looking at the two zoomed-in regions, one may further remark that the delineation of blood vessels given by the proposed method is sharper and more accurate.

In order to compare the performance of the two studied algorithms on the in-vivo data, two commonly employed quantitative measures are used. The first one is the resolution gain (RG) as suggested in [START_REF] Yu | A blind deconvolution approach to ultrasound imaging[END_REF]. It is defined as a ratio between the number of pixels of the normalized autocorrelations (values greater than 0.75, i.e. 3dB) of the ULM density images obtained respectively by the SVDbased algorithm (considered as a reference) and the proposed RPCA-based method. The calculated RG is 1.16, indicating a 16% improvement in the resolution of the proposed method over the SVD-based algorithm. The second one is contrast ratio (CR), a quantitative measure of contrast between the background and blood regions. A larger CR value indicates a better performance of ULM. For a pair of background and foreground patches, CR is defined as [START_REF] Rodriguez-Molares | The Generalized Contrast-to-Noise Ratio[END_REF]:

where µ R1 is the mean value of intensities in the blue rectangular patch R 1 which is kept fixed as a reference patch. µ R2 is the mean value of intensities in the green rectangular patch R 2 which is one of the non-overlapping patches of the image (see Fig. 1). Fig. 2 shows the boxplot of all the patchwise CR values of the two algorithms. One can see that the proposed method results in a much better contrast, producing a much higher median value (46.11 for the proposed PRCA-based ULM and 23.45 for the SVD-based ULM). It proves the consistency with the above qualitative observation about the ULM density image results of the two studied techniques.

V. CONCLUSIONS

In this paper, we introduced an efficient algorithm of ultrafast ultrasound localization microscopy (ULM) for hemodynamic maps of the brain microvasculature, based on improving the tissue filtering performance using the robust principal component analysis (RPCA)-based method. The proposed method enabled more accurate ULM results than the reference method. Assessments carried out on in vivo rat brain perfusion dataset showed the benefits of using the proposed algorithm. Future work could be devoted, for instance, to deep learning (DL) techniques for automatic hyperparameter tuning, alleviating a practical downside of current state-of-the-art methods.