Marie-Karelle Riviere

Jacques-Henri Jourdan

Sarah Zohar

dfcomb: an R-package for Phase I/II Trials of Drug Combinations

Keywords: Dose-finding, Phase I, Combination, Adaptive design, R package 1

Introduction

Most phase I cytotoxic dose-finding studies seek to establish a dose high enough to be able to observe potential efficacy while maintaining the toxicity rate within certain pre-defined acceptable limits. Phase I studies in oncology focus on determining the maximum tolerated dose (MTD) that will be used in further experiments [START_REF] Faries | Practical modifications of the continual reassessment method for phase I cancer clinical trials[END_REF][START_REF] Korn | A comparison of two phase I trial designs[END_REF][START_REF] O'quigley | Experimental designs for phase I and phase I/II dose-finding studies[END_REF]. However, nowadays, it is rare to find new agents that perform much better than existing therapeutic strategies. When combining two or more agents, there may be a synergistic effect in terms of efficacy. That is why investigators wish to increase overall anti-tumor action and survival by combining several agents: cytotoxics or targeted molecules or both. As a result, it is difficult to suppose that each agent will act independently in terms of toxicity.

When combining two agents, either a cytotoxic with a cytotoxic or a cytotoxic with a molecularly targeted agent (MTA), the underlying modeling assumptions diverge. Indeed, these two types of combinations differ in the shape of the dose-efficacy curves. For cytotoxic agents, a higher dose yields a greater response, but also leads to higher toxicity. However, the dose-efficacy relationship of the MTA may not follow the same monotonic pattern. The efficacy of some MTAs often increases at low dose levels and then plateaus (or approximately plateaus) at higher dose levels once a saturation level has been reached [START_REF] Druker | Perspectives on the development of a molecularly targeted agent[END_REF][START_REF] Gafter-Gvili | High-dose imatinib for newly di-agnosed chronic phase chronic myeloid leukemia patients-systematic review and metaanalysis[END_REF][START_REF] Ellis | Antiangiogenic therapy: more promise and, yet again, more questions[END_REF][START_REF] Morgan | Dynamic contrast-enhanced magnetic resonance imaging as a biomarker for the pharmacological response of PTK787/ZK 222584, an inhibitor of the vascular endothelial growth factor receptor tyrosine kinases, in patients with advanced colorectal cancer and liver metastases: results from two phase I studies[END_REF]. Recently, we proposed two novel statistical designs for each combination type [START_REF] Riviere | A Bayesian dose-finding design for drug combination clinical trials based on the logistic model[END_REF][START_REF] Riviere | A Bayesian dose finding design for clinical trials combining a cytotoxic agent with a molecularly targeted agent[END_REF]. In practice, these methods are based on Bayesian approaches that require high statistical and programming skills.

Many of the newly developed statistical methods for combination designs do not include software for their implementation. Some exceptions include POCRM [START_REF] Wages | pocrm: an R-package for phase I trials of combinations of agents[END_REF] and PIPE [START_REF] Mander | A product of independent beta probabilities dose escalation design for dual-agent phase I trials[END_REF]. As a result, many combination studies are still conducted with methods used for single agents [START_REF] Riviere | Designs of drugcombination phase I trials in oncology: a systematic review of the literature[END_REF].

In this article, we present a new R [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] package dfcomb, which provides functions for the conduct of a prospective single trial and simulation studies of the combination of two cytotoxics or a cytotoxic and an MTA in phase I dose-finding trials [START_REF] Riviere | A Bayesian dose-finding design for drug combination clinical trials based on the logistic model[END_REF][START_REF] Riviere | A Bayesian dose finding design for clinical trials combining a cytotoxic agent with a molecularly targeted agent[END_REF]. The primary features of dfcomb are as follows:

• Generates simulation replicates of phase I/II clinical trials for combination studies of two cytotoxic agents or a cytotoxic with an MTA.

• Determines the next or recommended combination in a phase I/II combination clinical trial.

The primary objective of this paper is to present the features of each method and to describe how to use the R dfcomb package in practice. Section 2 describes the models and computational methods of the package. Section 3 presents the R functions with simulations and examples of two dataset clinical trials, and in section 4 we conclude with a discussion.

Computational methods

Combination of two agents where toxicity and efficacy increase with both agents (CombIncrease_sim and CombIncrease_next)

Riviere et al [START_REF] Riviere | A Bayesian dose-finding design for drug combination clinical trials based on the logistic model[END_REF] proposed a Bayesian phase I dose-finding design for clinical trials combining two agents where the toxicity of both agents is assumed to increase with dose levels. As in classic phase I design, the design only accounted for toxicity because the efficacy was assumed to be correlated with toxicity and to increase with dose levels. Therefore, the aim was to determine the maximum tolerated dose (MTD) in the context of the two-drug combination, defined as the combination associated with a toxicity closest to a maximum pre-defined doselimiting toxicity (DLT) rate.

The combination toxicity probability at dose (j, k), π j,k , composed of dose level j of agent 1 (j = 1, . . . , J) and dose level k of agent 2 (k = 1, . . . , K), was modeled using a 4-parameter logistic model:

logit(π j,k ) = β 0 + β 1 u j + β 2 v k + β 3 u j v k ,
where u j and v k denote the standardized doses of agents 1 and 2, β 0 , β 1 , β 2 and β 3 are unknown parameters such that β 1 > 0, β 2 > 0, ∀k, β 1 + β 3 v k > 0 and ∀j, β 2 + β 3 u j > 0, ensuring that the toxicity probability increases with the dose levels of both agents, and intercept β 0 ∈ R. Gibbs sampling was used to sample in posterior distributions and estimate all posterior probabilities and parameters.

During the trial, determination of dose escalation and de-escalation to adjacent combinations was based on the uncertainty in P (π j,k < θ|data), that is on the probability that the current combination is safe enough. At the end of the trial, the MTD was selected as the tested combination with highest posterior probability, P (π j,k ∈ [θ -δ; θ + δ]), of being in a pre-defined interval around the targeted toxicity θ. We implemented an algorithm-based start-up phase to gather enough information at the beginning of the trial in order to have reliable estimates of posterior probabilities. Moreover, we added stopping rules in cases where all dose levels were too toxic or all dose levels were far from the targeted toxicity (and therefore potentially ineffective). The details of the method can be found in [START_REF] Riviere | A Bayesian dose-finding design for drug combination clinical trials based on the logistic model[END_REF].

Combination of two agents where the efficacy of one of the agents can plateau (CombPlateau_sim and CombPlateau_next)

Riviere et al [START_REF] Riviere | A Bayesian dose finding design for clinical trials combining a cytotoxic agent with a molecularly targeted agent[END_REF] proposed a Bayesian phase I/II dose-finding design for clinical trials combining two agents where the toxicity and efficacy of both agents are assumed to increase with dose levels, but the efficacy of one agent can plateau (e.g. molecularly targeted agent).

The toxicity is assumed to be a binary outcome (1 if DLT, 0 otherwise), and the toxicity probability at combination (j, k) was modeled using a logistic model:

logit(p jk ) = β 0 + β 1 u j + β 2 v k
where and u j and v k are standardized doses ascribed to the jth dose level of agent 1 and the kth dose level of agent 2, β 0 , β 1 , and β 2 are unknown parameters such that β 1 > 0 and β 2 > 0, so that toxicity monotonically increases with the dose levels of both agents 1 and 2.

In contrast, the efficacy response may require a longer follow-up time to be scored and was therefore considered as a time-to-event outcome. When combining an MTA with a cytotoxic agent, efficacy is expected to increase monotonically with the dose of the cytotoxic agent, but initially increases and then plateaus with the dose of the MTA, after the MTA reaches a level of saturation. We modeled the time to efficacy for the combination (j, k) using a proportional hazard model, λ jk (t), augmented with a plateau parameter τ :

λ jk (t) = λ 0 (t)exp{γ 1 w j + γ 2 (z k 1(k < τ ) + z τ 1(k ≥ τ ))},
where λ 0 (t) is the baseline hazard, 1(.) denotes the indicator function, and w j and z k are standardized doses. We assume that γ 1 > 0 and γ 2 > 0, and therefore efficacy increases monotonically with the dose of the cytotoxic agent 1. The plateau parameter τ is an integer between 1 and K and indicates at which dose level of agent 2 efficacy reaches a plateau. We assume an exponential distribution for the time to efficacy with a constant baseline hazard, i.e., λ 0 (t) = λ 0 , resulting in the following survival function for the time to efficacy:

S jk (t) = exp[-λ 0 t exp{γ 1 w j + γ 2 (z k 1(k < τ ) + z τ 1(k ≥ τ ))}].
Then, the response rate, q jk , at the end of the full follow-up time T for patients treated at the combination (j, k), is given by q jk = 1 -S jk (T ).

Gibbs sampling was used to sample in toxicity and efficacy posterior distributions and to estimate all posterior probabilities and parameters. The next optimal combination to be administered was selected as the admissible (safe: P (p jk > θ) < C T , and effective:

P (S jk (T ) > ξ) ≥ C E ;
where C T and C E are probability thresholds that need to be calibrated through a simulation study) combination with the highest estimate of efficacy, along with the lowest estimate of toxicity (if plateau reached).

We implemented an algorithm-based start-up phase to gather enough information at the beginning of the trial in order to have reliable estimates of posterior probabilities. If during the trial, all combinations were inadmissible, the trial was stopped without a recommendation of combination. The details of the method can be found in [START_REF] Riviere | A Bayesian dose finding design for clinical trials combining a cytotoxic agent with a molecularly targeted agent[END_REF]

R-functions

In this section, for each of the two methods described above, we present code examples and the corresponding outputs for: (1) simulation replicates of clinical trials, and (2) application to data of an ongoing clinical trial.

Combination of two agents where toxicity and efficacy increase with both agents (CombIncrease_sim and CombIncrease_next)

Simulation

We will illustrate the features of the CombIncrease_sim function. We considered 5 dose levels for the first agent (ndose_a1=5) and 3 dose levels for the second agent (ndose_a2=3).

Then we defined our scenario representing the true toxicity probabilities of the combination.

True toxicity probabilities should be entered in a matrix with agent 1 in rows and agent 2 in columns, with increasing toxicity probabilities with both row and column numbers, as for example:

p_tox <-matrix(c(0.05,0.10,0.15,0.30,0.45, 0.10,0.15,0.30,0.45,0.55, 0.15,0.30,0.45,0.50,0.60),nrow=5,ncol=3)

We chose a targeted toxicity probability (target) of 0.30, which is common in oncology due to the aggressiveness of the disease, and a targeted interval ([target_min;target_max]) of [0.20; 0.40]. We then provided initial guesses for toxicity probabilities for the dose levels of agent 1, and agent 2. Typically, if information from phase I clinical trials for these agents individually is available, it can be incorporated to set up the prior toxicity probabilities.

prior_tox_a1 <-c(0.12, 0.2, 0.3, 0.4, 0.5) prior_tox_a2 <-c(0.2, 0.3, 0.4)

In this example, if there is no early termination of the trial, a maximum total number of 20 cohorts (n_cohort) of 3 patients (cohort) should be included (that is 60 patients).

The toxicity outcome was binary (tite=FALSE), 1 if a DLT is observed, 0 otherwise. We performed 2000 replicates of phase I clinical trials (nsim).

In order to determine the combination of drugs for the next cohort of patients, probability thresholds c e and c d have to be defined. They are critical for the performance of the design as they control the dose escalation and de-escalation. When setting up a new clinical trial, the values of c e and c d should be carefully calibrated through simulation, to ensure good operating characteristics of the design. In practice, this can be performed as follows. First, define a set of representative dose-toxicity scenarios that may be encountered in the trial, and then perform simulations using different values of c e and c d to evaluate the performance of the design. This is a trial-and-error process and may involve repeatedly tuning the values of c e and c d based on the simulation results. The goal is to find the values of c e and c d that yield good overall performance across different scenarios (e.g., the percentage of correct selection of the MTD, the number of patients exposed to over-toxic combinations or undertoxic combinations). We fixed these settings to default values, which seems to give good performance in a large variety of scenarios [START_REF] Riviere | A Bayesian dose-finding design for drug combination clinical trials based on the logistic model[END_REF].

Finally, investigators may wish to stop a trial early when all combinations are estimated to be unacceptable in terms of toxicity. Therefore, we propose a stopping rule when the lowest combination is too toxic. In our example, to stop the trial early at least 6 patients should have been included in the lowest combination (n_min) and the probability that the estimated toxicity probability at that combination is above the target (0.30) should be greater than 0.95, meaning that the probability of over-dosing is high. Again, for a simpler use, these arguments were set by default to values giving, in general, good performances for inexperienced users.

Then users can choose to calibrate them by simulation studies. This results in the following call to the function: sim <-CombIncrease_sim(ndose_a1=5, ndose_a2=3, p_tox=p_tox, target=0.30, target_min=0.20, target_max=0.40, prior_tox_a1=prior_tox_a1, prior_tox_a2=prior_tox_a2, n_cohort=20, cohort=3, tite=FALSE, nsim=2000, c_e=0.85, c_d=0.45, c_stop=0.95, n_min=6, seed = 14061991)

A progress bar appears that indicates the progress of the number of simulations chosen.

The simulation can be long, and it is normal to wait some time before the bar changes.

Calculation time will clearly depend on the computer resources the user has at his or her disposal, but this example takes about 3 hours 45 minutes to run on a single portable computer with an Intel Core i7. If the call to the function is stopped before its end, the results for the current number of simulations ended are displayed. The generated output is as follows: Prior toxicity probabilities for agent 1:

[1] 0.12 0.20 0.30 0.40 0.50

Prior toxicity probabilities for agent 2: Toxicity is not a time-to-event but binary

Application

In a real clinical trial, the CombIncrease_next function determines the next combination to administer. As with CombIncrease_sim, the number of dose levels for each agent, targeted toxicities, initial guesses for toxicity probabilities, and minimum number of patients for early termination should be specified.

prior_tox_a1 <-c(0.12, 0.2, 0.3, 0.4, 0.5) prior_tox_a2 <-c(0.2, 0.3, 0.4)

In this example, we considered that the trial is still ongoing (final=FALSE), and currently 18 patients have been enrolled (pat_incl). For each of these 18 patients, two vectors corresponding to the dose levels of each agent administered to the patient should be given, as for example:

dose1 <-c(1,1,1,2,2,2,3,3,3,3,3,3,3,3,3,4,4,4) dose2 <-c(1,1,1,2,2,2,3,3,3,2,2,2,1,1,1,1,1,1)
In the following example, the last three patients have received dose 4 of agent 1 and dose 1 of agent 2. The toxicity outcome is binary (tite=FALSE), so the user have to input a vector of toxicity containing for each patient: 1 if a dose-limiting toxicity is observed and 0 otherwise.

toxicity <-c(0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,1)

If the toxicity is considered as a time-to-event (tite=TRUE), then two vectors containing the time-to-toxicity and the follow-up time for each patient as well as the value corresponding to the full follow-up time window should be given. The detailed example results in the following call to the function:

next_dose <-CombIncrease_next(ndose_a1=5, ndose_a2=3 , target=0.3, target_min=0.20, target_max=0.40, prior_tox_a1=prior_tox_a1, prior_tox_a2=prior_tox_a2, final=FALSE, pat_incl=18, dose_adm1=dose1, dose_adm2=dose2, tite=FALSE, toxicity=toxicity, n_min=6)
Another possible option include in_startup (set at TRUE by default). If FALSE the start-up phase is forced to be ended and thus estimation can be performed even if the clinical trial did not follow the start-up phase recommended by the original design. This code takes less than one second to run. The generated output is composed as the recalled number of

Simulation

We will illustrate the features of the CombPlateau_sim function. We considered 3 dose levels for the cytotoxic agent (ndose_a1) and 4 dose levels for the molecularly targeted agent (ndose_a2). We defined the true toxicity and efficacy probabilities of the combination. They should be entered in rows for the cytotoxic agent and in columns for the MTA, with increasing toxicity/efficacy probabilities with both row and column numbers:

p_tox <-t(matrix(c(0.10,0.15,0.30,0.45, 0.15,0.30,0.45,0.50, 0.30,0.45,0.55,0.65),nrow=4,ncol=3))

p_eff <-t(matrix(c(0.25,0.25,0.26,0.27, 0.40,0.41,0.41,0.42, 0.55,0.55,0.56,0.56),nrow=4,ncol=3))

We chose a maximum targeted toxicity probability of 0.30 (tox_max) and a minimum efficacy probability of 0.20 (eff_min). We then provided initial guesses for toxicity/efficacy probabilities for the dose levels of both agents:

prior_tox_a1 <-c(0.2, 0.3, 0.4) prior_eff_a1 <-c(0.3, 0.4, 0.5) prior_tox_a2 <-c(0.12, 0.2, 0.3, 0.4) prior_eff_a2 <-c(0.3, 0.4, 0.5, 0.59)

Again, if information from phase I clinical trials for these agents used as single-agents is available, it can be incorporated to set up the prior toxicity probabilities. In this example, if there is no early termination of the trial, a maximum sample size of 75 enrolled patients is needed (n), with a cohort of 3 patients throughout the trial (cohort_start, cohort). As efficacy is considered as a time-to-event outcome (time-to-progression), the full follow-up time to assess efficacy should be specified (time_full). In this example we chose a follow-up time of 7 weeks. The patient accrual is simulated using a Poisson process. We chose a rate of 0.28 (poisson_rate) meaning that on average 1 patient arrives every 3.5 weeks. Even if toxicity should be quickly ascertainable, it is possible to consider that before dose determination the cohort should be followed up for a minimum time to assess toxicity. Usually in oncology, 3 weeks are needed to evaluate toxicity. In the following example, we did not constrain a minimum follow-up time before dose determination and set cycle to 0. We performed 2000 replicates of phase I/II clinical trials (nsim).

In order to determine the combination of drugs for the next cohort of patients, toxicity and efficacy thresholds c tox and c ef f have to be defined. These probabilities are used to define combinations that are admissible in terms of toxicity and efficacy. A smaller c tox will be more stringent and require a safer combination in terms of toxicity, and a higher c ef f will require a more effective combination. The results show that the percentage of correct selection of the optimal combination (3, 1) is 73.25%. The trial was stopped at 2.95% because no combination was found that was acceptable in terms of both toxicity and efficacy.

Number of simulations: 2000

Cohort size for the start-up phase: 3

Cohort size for the model-based phase: 3

Total sample size: 75

Total patients accrued: 72.9

Toxicity upper bound: 0.3

Efficacy lower bound: 0.2

Prior toxicity probabilities for agent 1:

[1] 0.2 0.3 0.4

Prior toxicity probabilities for agent 2:

[1] 0.12 0.20 0.30 0.40

Prior efficacy probabilities for agent 1:

[1] 0.3 0.4 0.5

Prior efficacy probabilities for agent 2:

[1] 0.30 0.40 0.50 0.59 Toxicity threshold: 0.85 Efficacy threshold: 0.1

Toxicity is binary

Efficacy is a time-to-event Full follow-up time: 7

Minimum waiting time between two dose cohorts is of one toxicity cycle of: 0 toxicity <-c(0,0,0,0,0,1,0,0,0,0,0,1,0,0,1,0,0,0)

The full follow-up time for efficacy assessment was set to 7 weeks. As efficacy is assessed as a time-to-progression, for each patient, the follow-up time (maximum equal to the full follow-up time) as well as the time-to-progression are required. If no progression is observed for the patient, then the value should be set to +Inf. agents), or the efficacy can plateau with the dose of one agent (e.g. combination of a cytotoxic agent with a molecularly targeted agent). For each method, the package implements a single prospective clinical trial or simulation studies. In this package, several flexible inputs can be used, such as: (1) prior information on each agent gathered before the trial, (2) the cohort size of patients allocated at each combination level as well as the total sample size, (3) the targeted toxicity, (4) binary or censured data among other possible entries. In the last decade some innovative methods have been proposed for combination studies. However, the implementation of such complex approaches is difficult as most commercial software programs do not provide functions or macros for these approaches. Indeed, many early phase clinical trial methods are not used in practice as programs or scripts are not available. Our R package enables statisticians and physicians to implement adaptive model-based dose-finding designs when two anti-cancer agents are studied.

  True toxicity probabilities entered are recalled, and then the results show at each combination (1) the percentage of selection, (2) the mean number of patients included, and (3) the mean number of observed toxicities. Therefore, the percentage of correct selection of a true MTD (combination associated with the toxicity target of 30%) is 28.30 + 40.65 + 6.90 = 75.85%. The mean number of patients and toxicities of each combination can be used to evaluate the number of patients exposed to over-dosing or under-dosing.

3 . 2

 32 a time-to-event but binary Combination of two agents where the efficacy of one of the agents can plateau (CombPlateau_sim and CombPlateau_next)

3 )

 3 time_prog <-c(1.6,4.2,3.5,5.1,2.4,4.8,2.8,4.4,+Inf,3.9,+Inf,4.6,1.8,+Inf, 0.5, 5.4,2.8,+Inf) time_follow <-c(rep(7,15), 4.9, 3.1, 1.The example detailed above results in the following call to the function: next_dose <-CombPlateau_next(ndose_a1=3, ndose_a2=4, tox_max=0.30, eff_min=0.20, prior_tox_a1, prior_tox_a2, prior_eff_a1, prior_eff_a2, stage=0, in_startup=FALSE, cohort=3, pat_incl=18, dose_adm1=dose_adm1, dose_adm2=dose_adm2, toxicity=toxicity, time_full=7, time_prog=time_prog, time_follow=time_follow)This code takes two seconds to run. The output reports the number of patients included, observed toxicities and progressions at each combination according to the entries:

  When setting up a new clinical trial, the values of c tox and c ef f

	1 0.98 2.83 19.28
	Percentage of inconclusive trials:	2.95
	should be calibrated through simulation depending on the desired operating characteristics
	of the design. The following call to the function was performed:
	sim <-CombPlateau_sim(ndose_a1=3, ndose_a2=4, p_tox=p_tox, p_eff=p_eff,
	tox_max=0.30, eff_min=0.20, prior_tox_a1=prior_tox_a1,
	prior_tox_a2=prior_tox_a2, prior_eff_a1=prior_eff_a1,
	prior_eff_a2=prior_eff_a2, n=75, cohort_start=3, cohort=3, time_full=7,

The trials did not stop early and included the total sample size of 60 patients. Other entries are also recalled as follows:

Percentage of inconclusive trials: 0

The minimum number of cohorts to stop the trial is: The estimated posterior probabilities for toxicity, under/over-dosing are reported. Prior toxicity probabilities for agent 1:

[1] 0.12 0.20 0.30 0.40 0.50

Prior toxicity probabilities for agent 2: poisson_rate=0.28, cycle=0, nsim=2000, c_tox=0.85, c_eff=0.10, seed = 2174892, threads=0)

The seed for simulations can be chosen. Finally, as these simulations are computationally intensive, the run time can be long. The argument threads is used to specify the number of threads of the computer to use; by default 0 uses all threads. As this design deals with both toxicity and efficacy, simulations in general take longer to execute than for the previous design. We run this code in a little over 2 days (52 hours).

The generated output is first composed of the recalled toxicity and efficacy probabilities:

True toxicities: The mean trial duration was 268 weeks, that is about 5 years, if only one center is recruiting, with on average one patient enrolled every 3.5 weeks. Entries and accumulated data are reported.

Application

When used for a real clinical trial, the CombPlateau_next function determines the next combination to administer. As with CombPlateau_sim, the number of dose levels for each agent, maximum toxicity, minimum efficacy, and initial guesses for toxicity/efficacy probabilities should be specified: stage is an integer with value 0 if less than half of the total sample size has been included, 1 if more than half of the total sample size has been included but the trial is still ongoing, and 2 if the trial is over and the dose recommendation should be made. We considered that the trial is ongoing with less than half of the total sample size included, setting stage to 0, and the start-up phase is finished (in_startup=FALSE). After the start-up phase, patients were included by cohorts of 3 patients (cohort). Currently, 18 patients have been enrolled (pat_incl). For each patient, two vectors corresponding to the dose levels of each agent administered to the patient should be given:

A vector containing the binary toxicity outcome for each patient has to be defined. Then estimated toxicities and efficacies leading to recommendation of combinations are given:

Toxicity prob: Prior toxicity probabilities for agent 1:

[1] 0.2 0.3 0.4

Prior toxicity probabilities for agent 2:

[1] 0.12 0.20 0.30 0.40

Prior efficacy probabilities for agent 1:

[1] 0.3 0.4 0.5

Prior efficacy probabilities for agent 2: Minimum waiting time between two dose cohorts is of one toxicity cycle of: 0

Conclusions

The dfcomb package implements two methods for phase I clinical trials in oncology for the combination of two agents depending on the underlying assumptions: either both toxicity and efficacy increase with the dose of each agent involved (e.g. combination of two cytotoxic


K Rojkova

E Volle

M Urbanski

F Humbert

F Dell'acqua

Michel Thiebaut De Schotten

Acqua F Dell'

Thiebaut De Schotten

Atlasing the frontal lobe connections and their variability due to age and education: a spherical deconvolution tractography study

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.

The study of frontal lobe connections is essential to understanding of the role and the organization of the frontal lobes' distinct subregions. The functions of the frontal subregions emerge from an interaction with other cortical areas, on which it exerts its control, and from which it receives information to integrate. These interactions are supported by short and longrange white matter connections. While most of the recent data and theories regarding the functional organization of the frontal lobes have been supported by functional imaging, its anatomical connectivity has yet only been scarcely investigated.

In 1985, a breakthrough came with the emergence of a new magnetic resonance imaging (MRI) modality: the diffusion-weighted MRI (Le [START_REF] Bihan | Imagerie de diffusion in-vivo par résonance magnétique nucléaire[END_REF]. This non-invasive technique offers a unique approach to investigate in-vivo the structure of brain tissues by measuring the diffusion of water molecules along different directions. This exciting technical advance not only reproduced reliably known anatomy in the living human brain [START_REF] Catani | Virtual in vivo interactive dissection of white matter fasciculi in the human brain[END_REF]Catani and Thiebaut de Schotten 2008), but also refined the description of large fiber bundles [START_REF] Catani | Perisylvian language networks of the human brain[END_REF]. Moreover, it also gave raise to studies of intersubject variability [START_REF] Ciccarelli | Diffusion tractography based group mapping of major white-matter pathways in the human brain[END_REF]Thiebaut de Schotten et al. 2011;[START_REF] Verhoeven | Construction of a stereotaxic DTI atlas with full diffusion tensor information for studying white matter maturation from childhood to adolescence using tractography-based segmentations[END_REF][START_REF] Wakana | Reproducibility of quantitative tractography methods applied to cerebral white matter[END_REF], asymmetries [START_REF] Barrick | White matter pathway asymmetry underlies functional lateralization[END_REF]Catani et al. 2007), behavioral correlates [START_REF] Barrick | White matter pathway asymmetry underlies functional lateralization[END_REF]Catani et al. 2007;[START_REF] Johansen-Berg | Integrity of white matter in the corpus callosum correlates with bimanual co-ordination skills[END_REF]; Thiebaut de Schotten et al. 2014a) and pathological differences [START_REF] Craig | Altered connections on the road to psychopathy[END_REF]Thiebaut de Schotten et al. 2014a). While several fiber populations cross almost everywhere in the brain (Dell'Acqua et al. 2013;[START_REF] Jeurissen | Investigating the prevalence of complex fiber configurations in white matter tissue with diffusion magnetic resonance imaging[END_REF]) standard diffusion tensor tractography models follow the principal direction of fibers, leading to incomplete or erroneous reconstructions of white matter features (Thiebaut de Schotten et al. 2011) that may bias tract-specific measurements [START_REF] Dell'acqua | Structural human brain networks: hot topics in diffusion tractography[END_REF][START_REF] Vos | The influence of complex white matter architecture on the mean diffusivity in diffusion tensor MRI of the human brain[END_REF]. New developments in tractography and diffusion modeling can circumvent this issue by extracting the orientation distribution of different populations of fibers within the same voxel (Dell'acqua et al. 2010;[START_REF] Descoteaux | Regularized, fast, and robust analytical Q-ball imaging[END_REF]Tournier et al. 2004;[START_REF] Tuch | Diffusion MRI of complex neural architecture[END_REF][START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF]. These new methods allow for the depiction of tracts that better correspond to postmortem dissections [START_REF] Dauguet | Comparison of fiber tracts derived from in-vivo DTI tractography with 3D histological neural tract[END_REF][START_REF] Dauguet | 3D histological reconstruction of fiber tracts and direct comparison with diffusion tensor MRI tractography[END_REF][START_REF] Lawes | Atlasbased segmentation of white matter tracts of the human brain using diffusion tensor tractography and comparison with classical dissection[END_REF]Thiebaut de Schotten et al. 2011) and simian axonal tracing (Thiebaut de Schotten et al. 2012) with good reliability (Kristo et al. 2013b). More recently, a novel nomenclature of the frontal lobe U-shaped connection was made possible through preliminary results that combined post-mortem Klingler dissection [START_REF] Klingler | Erleichterung der makroskopischen Präparation des Gehirn durch den Gefrierprozess[END_REF] with high resolution spherical deconvolution tractography (Catani et al. 2012). However, intersubject variability and aging effects could not be assessed due to the small sample size of participants studied. New tract-specific measurements taking into account the crossing of fibers also followed these advances in diffusion imaging tractography (e.g. Apparent Fiber Density or Hindrance Modulated Orientational Anisotropy (HMOA) [START_REF] Raffelt | Apparent Fibre Density: a novel measure for the analysis of diffusion-weighted magnetic resonance images[END_REF]Dell'Acqua et al. 2013). These new measurements are specific to the direction followed by the tractography; therefore providing a true-tract specific index to characterize white matter diffusion when crossing. These indices could offer a more precise characterization of the underlying microstructural organization, and consequently, be more sensitive to age-related differences of individual tracts when compared to traditional diffusion tensor imaging indices, such as fractional anisotropy (FA).

The study of the frontal connections is also of particular interest for the neurosciences of aging. With aging, neuronal loss and small vessel alteration lead to progressive and subtle white matter changes associated with cognitive decline in the elderly [START_REF] Pantoni | Cerebral small vessel disease: from pathogenesis and clinical characteristics to therapeutic challenges[END_REF]Xiong and Mok 2011). Cognitive decline affects predominantly many executive functions but not all [START_REF] Geerligs | Brain mechanisms underlying the effects of aging on different aspects of selective attention[END_REF]) and brain changes seem to distribute unevenly, affecting predominantly the frontal region [START_REF] Bishop | Neural mechanisms of ageing and cognitive decline[END_REF][START_REF] Curiati | Brain structural variability due to aging and gender in cognitively healthy Elders: results from the Sao Paulo Ageing and Health study[END_REF][START_REF] Draganski | Regional specificity of MRI contrast parameter changes in normal ageing revealed by voxel-based quantification (VBQ)[END_REF][START_REF] Giorgio | Age-related changes in grey and white matter structure throughout adulthood[END_REF][START_REF] Good | Cerebral asymmetry and the effects of sex and handedness on brain structure: a voxel-based morphometric analysis of 465 normal adult human brains[END_REF][START_REF] Moscovitch | Memory and Working-with-Memory: A Component Process Model Based on Modules and Central Systems[END_REF][START_REF] Raz | Aging, sexual dimorphism, and hemispheric asymmetry of the cerebral cortex: replicability of regional differences in volume[END_REF][START_REF] Raz | Neuroanatomical and cognitive correlates of adult age differences in acquisition of a perceptual-motor skill[END_REF]West 1996) suggesting that latest developed brain areas are especially susceptible to show an accelerated aging [START_REF] Tamnes | Brain Maturation in Adolescence and Young Adulthood: Regional Age-Related Changes in Cortical Thickness and White Matter Volume and Microstructure[END_REF]. Previous cross sectional studies using tract-specific measurements confirmed this "frontal lobe hypothesis" for the white matter by revealing a slow decrease with aging in the FA for frontal callosal tracts [START_REF] Hsu | Gender differences and age-related white matter changes of the human brain: a diffusion tensor imaging study[END_REF]Hasan et al. 2009b;[START_REF] Lebel | Age-related regional variations of the corpus callosum identified by diffusion tensor tractography[END_REF][START_REF] Bastin | Tract shape modelling provides evidence of topological change in corpus callosum genu during normal ageing[END_REF][START_REF] Michielse | Selective effects of aging on brain white matter microstructure: a diffusion tensor imaging tractography study[END_REF]) and for long tracts connecting the frontal lobe [START_REF] Jones | Age effects on diffusion tensor magnetic resonance imaging tractography measures of frontal cortex connections in schizophrenia[END_REF]Hasan et al. 2009a;[START_REF] Hasan | Quantification of the spatiotemporal microstructural organization of the human brain association, projection and commissural pathways across the lifespan using diffusion tensor tractography[END_REF], affecting more prominently the frontal portion of these tracts [START_REF] Davis | Assessing the effects of age on long white matter tracts using diffusion tensor tractography[END_REF]). The same result was corroborated with voxel-based statistics showing a negative correlation between aging and FA in specific but not all areas of the frontal lobes [START_REF] Giorgio | Age-related changes in grey and white matter structure throughout adulthood[END_REF][START_REF] Madden | Adult age differences in the functional neuroanatomy of visual attention: a combined fMRI and DTI study[END_REF][START_REF] Madden | Diffusion tensor imaging of adult age differences in cerebral white matter: relation to response time[END_REF][START_REF] O'sullivan | Evidence for cortical "disconnection" as a mechanism of age-related cognitive decline[END_REF][START_REF] Phillips | Superficial white matter: effects of age, sex, and hemisphere[END_REF][START_REF] Raz | The aging brain observed in vivo: differential changes and their modifiers, in: Cognitive Neuroscience of Aging: Linking Cognitive and Cerebral Aging[END_REF][START_REF] Salat | Prefrontal gray and white matter volumes in healthy aging and Alzheimer disease[END_REF]). Nevertheless, it is important to note that the "frontal lobe hypothesis" has also been contested by a two years longitudinal study showing no evidence of an accelerated decline in the frontal lobe region in the aging population [START_REF] Barrick | White matter structural decline in normal ageing: a prospective longitudinal study using tract-based spatial statistics[END_REF]. Aging may affects some, but not all, tracts of frontal white matter. This explains the contrasted results reported in the literature. Age-related differences are not taken into account in atlases built from populations with a restricted age-range; this bias may hamper the use of these atlases in various clinical populations. One way to circumvent this issue is to build a normative global atlas based on a healthy population with a large age range. The effects of age may also differ according to some factors, such as compensatory networks [START_REF] Stern | Brain networks associated with cognitive reserve in healthy young and old adults[END_REF][START_REF] Stern | A common neural network for cognitive reserve in verbal and object working memory in young but not old[END_REF] or increased functional connectivity in existing networks [START_REF] Bastin | Cognitive reserve impacts on inter-individual variability in resting-state cerebral metabolism in normal aging[END_REF], providing a partial protection called "cognitive reserve" [START_REF] Stern | What is cognitive reserve? Theory and research application of the reserve concept[END_REF][START_REF] Stern | Cognitive reserve[END_REF]. Higher education seems to also be an important factor delaying aging affects on the brain and cognitive decline [START_REF] Brayne | Education, the brain and dementia: neuroprotection or compensation?[END_REF][START_REF] Coffey | Relation of education to brain size in normal aging: implications for the reserve hypothesis[END_REF]) and preliminary evidences suggest structural connectivity changes associated with cognitive reserve in elderly [START_REF] Fischer | Association of structural global brain network properties with intelligence in normal aging[END_REF]). However it is unknown whether specific tracts are supporting these changes.

Therefore, the aim of our study was to map, in a sample of healthy participants, all the association, projection and short U-shaped tracts previously described in the frontal lobe. We also estimated the effect of age and level of education on these tracts using measures of the volume, FA and HMOA. The final atlas of the human frontal lobe can be used for future clinicalneuroanatomical correlation in patients of various age ranges and levels of education.

MATERIAL AND METHODS

Participants

This study was approved by the local Ethics Committee. Participants were recruited via advertising in the Salpêtriere hospital and on the web site www.risc.cnrs.fr. All participants were right-handed healthy adults with no history of neurological or psychiatric disorder, no cognitive complaints, and no cognitive impairments or depression, as assessed using translated versions of the Mini Mental state [START_REF] Folstein | Mini-mental state". A practical method for grading the cognitive state of patients for the clinician[END_REF]) and the MADRS [START_REF] Montgomery | A new depression scale designed to be sensitive to change[END_REF].

Participants were not included if their MMSE score was lower than 27 and if they have history of neurological or psychiatric symptoms.

All subjects had an MRI scan including 3D T1-weighted images in addition to diffusion images (acquisition, preprocessing, rules for dissections as well as the approach employed to map white matter tracts are reported in supplementary material). All brain images were examined by a neuroradiologist. Millimetric T1-weighted and diffusion weighted images did not show any significant signal abnormalities evocative of a small vessel disease or of an evolving neurological disease. Subjects with abnormalities on MRI were excluded.

Informed and written consent to participate in this study was provided by 57 right-handed volunteers. Participants were excluded from the analysis for medical reasons (anomalies on neuropsychological testing or on the brain MRI) or for head movements during the MRI acquisition. The average age of the 47 remaining participants (24 males and 23 females) was 45.45 years (±14.79 years; aged between 22 and 71 years) and their education level was 15.36 years of education (±3.00 years; range between 10 and 26 years), noting that a French bachelor's degree corresponds to 12 years of education. Subjects from our sample had relatively high education levels, independently of their gender (t (45) = 1.203; p = 0.235) but with a trend of correlation with age, younger participants having a slightly but not significantly higher education level than older participants (r= -0.259; p = 0.079). Therefore partial correlation statistics were carried on controlling for age or education when appropriate.

Statistical analysis

In order to test whether white matter differences distribute unevenly, affecting predominantly the frontal region, we used SPSS software (SPSS, Inc., Chicago, IL, United States of America) to carry on a repeated ANOVA between tracts emerging from the frontal lobe and tracts not involving the frontal lobe using age and education as covariates.

For each of the following statistics, we employed a False Discovery Rate (FDR) correction [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] for multiple comparisons; this is available as a tool on the Signed Differential Mapping website (http://www.sdmproject.com/).

Gaussian distribution of the data using the Shapiro-Wilk test [START_REF] Shapiro | An Analysis of Variance Test for Normality (Complete Samples)[END_REF] was not confirmed for all variables in our group of participants. Therefore, two-tailed partial Spearman rank correlation coefficients [START_REF] Spearman | The proof and measurement of association between two things[END_REF] were performed between the measurements of each dissected single tract and age, controlling for the education level of the participants. Similarly, partial Spearman ranking correlation coefficients were performed between the measurements of each dissected single tract and the education level controlling for age (FDR corrected for 55 tracts).

RESULTS

We successfully identified and mapped, in the MNI reference space, 55 frontal tracts: 16 interlobar association tracts, one commissural tract, eight frontal projection tracts and 30 frontal short U-shaped tracts for both hemispheres (Fig. 1). Maps are provided as an online digital supplementary material (https://www.dropbox.com/s/g0vlklnk60up73y/Atlasmaps.zip) and projections of these tracts are provided in Fig. 2, 3 and4. Percentage of success for the reconstruction of each tract and tract-specific measurements are tabulated in supplementary table 1. This percentage was calculated using as follows:

Percentage of success = Number of successful reconstructions of the tract/number of participants X 100.

Frontal interlobar association tracts

Superior longitudinal fasciculus (SLF)

We reconstructed a system, composing of three parallel longitudinal branches passing through the dorsolateral portion of the white matter, situated above the ventricles and the lateral sulcus, and connecting the parietal and the frontal lobes. The superior branch (SLF I; n°1 in Fig. 1 and Fig. 2) runs from the superior parietal lobule and the precuneus (Brodman Area, BA7 and BA5) to the superior frontal gyrus (BA6, BA8, BA9 up to BA10). More ventrally, the middle branch (SLF II; n°2 in Fig. 1 and Fig. 2) connects the angular gyrus (BA39) to the posterior regions of the middle frontal gyrus (BA6). Few projections continue further into the middle frontal gyrus up to BA46. The inferior branch (SLF III; n°3 in Fig. 1 and Fig. 2) originates in the supramarginal gyrus (BA40) and terminates in the pars opercularis (BA44), triangularis (BA45) and the inferior frontal gyrus (BA47). Fig. 1 andFig. 2) We identified a long medial fiber bundle running within the cingulate gyrus, all around the corpus callosum, that contains fibers of different lengths. The longest fibers run from the parahippocampal gyrus and the uncus (BA34) to the medial portion of the orbitofrontal (BA11) cortex. The cingulum also contains short U-shaped fibers that connect the medial frontal (BA32, BA8, BA9, BA10) and parietal (BA31, BA5, BA7) to different portions of the cingulate gyrus.

Cingulum (n°4 in

Uncinate (n°5 in Fig. 1 and Fig. 2)

The fibers of the uncinate originate from the temporal pole (BA38), parahippocampal gyrus, uncus (BA34) and amygdala. After a U-turn, they enter the external capsule then reach the lateral and ventral orbital cortex (BA11, BA12) and frontal pole (BA10).

Arcuate fasciculus

The fronto-temporal portion of the arcuate fasciculus connects the pars opercularis (i.e.

Broca area, BA44) to the auditory cortex (BA22) and the posterior portion of the middle and inferior temporal gyri (BA21, BA37). It forms the long segment (LS; n°6 in Fig. 1 and Fig. 2) of the arcuate fasciculus, arching around the sylvian fissure. Shorter fronto-parietal fibers link the ventral portion of the precentral gyrus (BA6) to the postcentral gyrus (BA1) and the supramarginal gyrus (BA40). These fibers form the anterior segment (AS; n°7 in Fig. 1 and Fig.

2) of the arcuate fasciculus.

Inferior fronto-occipital fasciculus (IFOF, n°8 in Fig. 1 and Fig. 2)

The inferior fronto-occipital fasciculus is a long-ranged bow-tie-shaped tract connecting the occipital lobe to the frontal lobe. It originates from the occipital pole (BA18), the lingual gyrus and the inferior and middle occipital gyri (BA19). As it leaves the occipital lobe, IFOF narrows and its fibers gather at the level of the extreme capsule. As it enters the frontal lobe, its fibers spread to the inferior frontal gyrus. Its ventral fibers reach the orbital gyrus (BA11) and the frontal pole (BA10), whereas its dorsal fibers terminate in the anterior superior frontal gyrus (most anterior part of BA9 and dorsal BA10).

Frontal commissural and projection tracts

Commissural and projection tracts are difficult to follow with tractography as they intertwine with one another and cross association fibers. Despite spherical deconvolution being designed to circumvent this fiber-crossing issue, their complex configuration may occasionally lead to an underestimation of their true extent. Note that the validity and limitations of these results will be highlighted in the discussion.

Frontal Corpus callosum (n°9 in Fig. 1 and Fig. 3)

The frontal portion of the corpus callosum occupies the rostrum, the genu and most of the body of the corpus callosum and projects on to the whole surface of the frontal lobe except for the pars orbitalis (BA 47). These projections also show a reduced probability for the pars opercularis (BA 44) and the frontal eye field (BA 6).

Corticospinal tract (CST, n°10 in Fig. 1 and Fig. 3)

The corticospinal tract consists of a fan-shaped set of fibers, passing through the posterior limb of the internal capsule and projecting from the precentral gyrus (BA4) to the spinal cord. We note that in a third of participants we obtained its projection onto the whole homunculus in the precentral gyrus; for the remainder of the participants tractography confined its projections to the hand and the lower limb regions.

Fronto-thalamic (anterior thalamic radiations) and fronto-striatal projections (n°11 and 12 in

Fig. 1 and Fig. 3)

Fronto-thalamic and fronto-striatal projections pass through the anterior limb of the internal capsule and projects on the whole surface of the frontal lobe. These projections also show a reduced probability for the ventral portion of the precentral gyrus and the middle frontal gyrus (BA 6 and 46).

Fronto-pontine projections (n°13 in Fig. 1 and Fig. 3) Fronto-pontine fibers pass through the genu and the anterior limb of the internal capsule and projects on to the whole surface of the frontal lobe. As for the fronto-thalamic and frontostriatal projections, we note a reduced probability of the fronto-pontine projection in the precentral gyrus and the middle frontal gyrus (BA 6 and 46).

Frontal U-shaped tracts

Fronto-parietal U-tracts

We mapped a set of five tracts running beneath the central sulcus from the parietal postcentral gyrus (BA1) to frontal precentral gyrus (BA4), which are organized from dorsal to ventral as follows: paracentral U-shaped tract (n°14 in Fig. 1 and Fig. 4), ventral face U-shaped tract (n°18 in Fig. 1 and Fig. 4), and hand superior, middle and inferior U-shaped tracts (n°15 to 17 in Fig. 1 and Fig. 4).

Frontal aslant tract (FAT) (n°19 in Fig. 1 and Fig. 3)

The frontal aslant tract is a large tilted tract linking the supplementary and presupplementary motor area (BA6) with the frontal operculum (BA44) and the posterior portion of the pars triangularis (BA45).

Fronto-insular tracts (FITs) (n°20 to 24 in Fig. 1 and Fig. 3)

A system composing of five U-shaped tracts dives from the frontal lobe into the external capsule to reach the insula. The four most anterior tracts, 1 to 4, link respectively the pars orbitalis or BA47 (FIT 1; n°20 in Fig. 1 and Fig. 4), the pars triangularis or BA45 (FIT 2; n°21 in Fig. 1 and Fig. 4), the pars opercularis or BA44 (FIT 3; n°22 in Fig. 1 and Fig. 4) and the precentral gyrus or BA6 (FIT 4; n°23 in Fig. 1 and Fig. 4) of the frontal lobe with the anterior insula (i.e. anterior short gyrus and middle short gyrus, see [START_REF] Cerliani | Probabilistic tractography recovers a rostrocaudal trajectory of connectivity variability in the human insular cortex[END_REF] for a detailed description of the anatomy of the insula). The most posterior tract (FIT 5; n°24 in Fig. 1 and Fig. 4) connects subcentral gyrus (BA43) to the posterior insula (i.e. anterior long gyrus).

Frontal longitudinal system

Two parallel chains of U-shaped tracts run longitudinally in the frontal lobe: the frontal superior longitudinal tract (FSL; n°25 in Fig. 1 and Fig. 4) and frontal inferior longitudinal tract (FIL; n°26 in Fig. 1 and Fig. 4).

The longest fibers of frontal superior longitudinal tract connect the precentral gyrus to the anterior portion of the superior frontal gyrus, while its shortest fibers link the middle frontal gyrus with the superior frontal gyrus through a series of a chain-like U-shaped fibers.

The longest fibers of frontal inferior longitudinal tract connect the precentral gyrus to the anterior portion of the middle frontal gyrus, while its shortest fibers link the dorsal portion of the pars opercularis and triangularis with the middle frontal gyrus through a series of chain-like Ushaped fibers.

Fronto-orbito-polar and fronto-marginal tracts

Two intralobar tracts can be observed within the frontal pole. The fronto-orbito-polar tract (FOP; n°27 in Fig. 1 and Fig. 4) connects the posterior orbital gyrus (BA12) to the anterior orbital gyrus (BA11) and the ventromedial region of the frontal pole.

The fronto-marginal tract (FMT; n°28 in Fig. 1 and Fig. 4) runs beneath the frontomarginal sulcus from the medial to the lateral region of the frontal pole (BA10).

Age and level of education effect

Repeated ANOVA between measurement of streamlines emerging from the frontal lobe and streamlines not involving the frontal lobe using age and education as covariates revealed a significant interaction between HMOA measures of frontal and non-frontal tracts and the age of the participants (F (1,44) =6.968; p<0.011). This result demonstrates that microstructural differences related to age are stronger for the tracts emerging from the frontal lobe compared to other tracts. This effect was not significant for FA measure or relative volume measures.

Results presented in the following section are false discovery rate (FDR) corrected for multiple comparisons. Uncorrected p values are reported in supplementary table 2. Tract specific volume measurements did not correlate significantly with aging.

Partial correlation analyses revealed that aging was significantly associated with a decrease of FA (Spearman rank coefficient partial correlation δ = -0.519, p < 0.001) and HMOA (δ = -0.553, p < 0.001) in the frontal projections of the corpus callosum. HMOA measures also decreased with aging in the left hemisphere for the frontal inferior longitudinal fasciculus (δ = -0.563, p < 0.001), the frontal orbito-polar tract (δ = -0.552, p < 0.001), fronto-insular tract 5 (δ = -0.379, p = 0.047), the corticospinal tract (δ = -0.442, p < 0.014) and fronto-thalamic projections (δ = -0.499, p = 0.004) (Fig. 5 and6). No other partial correlation with age was reported as significant.

Aging was also associated with a decrease of HMOA in the right frontal lobe including the SLF I (δ = -0.452, p = 0.015), SLF III (δ = -0.610, p < 0.001) branches of the superior longitudinal fasciculus, the corticospinal tract (δ = -0.380, p < 0.046) and fronto-thalamic projections (δ = -0.518, p = 0.003) (Fig. 5 and6).

Tract specific volume, FA and HMOA correlations with education level did not survive FDR correction for multiple comparisons.

DISCUSSION

In this study we employed advanced spherical deconvolution tractography to dissect 55 white matter pathways within the frontal lobes. We produced a probabilistic atlas of the interlobar, commissural, projection and short U-shaped tracts in a population with a large age range. All tracts showed a high inter-subject spatial reproducibility (downloadable as supplementary material; https://www.dropbox.com/s/g0vlklnk60up73y/Atlasmaps.zip). We also replicated correlations between age and microstructural measures of the corpus callosum, bilateral corticospinal tract and fronto-thalamic projections. For the first time, the use of spherical deconvolution combined with the use of tract specific HMOA allowed us to report correlations between age and the frontal inferior longitudinal fasciculus, the fronto-insular tract 5 and orbitopolar tract in the left hemisphere and the SLF I and III in the right hemisphere. Education did not correlate with any of the tract-specific measurements. While our maps can be of interest for clinical neuroanatomical correlations and functional imaging studies, the correlation with age provide a normative measure for future studies on diseases affecting various age ranges.

Probabilistic atlas of the interlobar, commissural, projection and short U-shaped tracts

Our frontal lobe atlas is a probabilistic replication in the stereotaxic space (MNI152) of the original anatomical post-mortem and tractography dissection performed on a single brain and reported by Catani et al. (2012). Compared to previous work [START_REF] Mori | Stereotaxic white matter atlas based on diffusion tensor imaging in an ICBM template[END_REF][START_REF] O'donnell | Automatic tractography segmentation using a high-dimensional white matter atlas[END_REF][START_REF] Lawes | Atlasbased segmentation of white matter tracts of the human brain using diffusion tensor tractography and comparison with classical dissection[END_REF]Thiebaut de Schotten et al. 2011), this atlas provides probabilistic maps of the three branches of the superior longitudinal fasciculus and many intralobar and U-shaped connections, as well as a full reconstruction of most of the crossing fibers. These maps can be of interest to identify tracts connecting functional activations (Thiebaut It may also provide bases for future studies aiming to clarify the very debated anatomo-functional organization of the prefrontal cortex for executive functions and cognitive control.

Classically, patient studies demonstrate that the prefrontal cortex is critical to executive functions, and adaptive goal-directed behaviors [START_REF] Azuar | Testing the model of caudo-rostral organization of cognitive control in the human with frontal lesions[END_REF][START_REF] Burgess | Mesulam's frontal lobe mystery re-examined[END_REF][START_REF] Stuss | Principles of Frontal Lobe Function[END_REF][START_REF] Volle | The functional architecture of the left posterior and lateral prefrontal cortex in humans[END_REF]. Influent cognitive models such as the working memory model postulate a specialization of retro-rolandic and of prefrontal regions along the superior-inferior axis. In the posterior prefrontal regions, dissociation between the processing of spatial (in the superior part) and verbal (in the inferior part) information has been reported in both functional MRI [START_REF] Curtis | Success and failure suppressing reflexive behavior[END_REF][START_REF] Sakai | Prefrontal interactions reflect future task operations[END_REF][START_REF] Sala | Functional topography of a distributed neural system for spatial and nonspatial information maintenance in working memory[END_REF][START_REF] Smith | Storage and executive processes in the frontal lobes[END_REF], and lesion studies [START_REF] Volle | The functional architecture of the left posterior and lateral prefrontal cortex in humans[END_REF][START_REF] Volle | A new era for lesion-behavior mapping of prefrontal functions[END_REF]. The superior and inferior longitudinal U-shaped tracts evidenced in the current study may support this spatial and verbal specialization respectively within the frontal lobe. Posteriorly, separate parietal and temporal inputs to these systems via the superior longitudinal and the arcuate fasciculi may also support this functional segregation. The examination of direct relationship between anatomical tracts and functional dissociations would be of high interest for future studies in order to disambiguate the organization of prefrontal cortex for executive functions.

The longitudinal U-shaped tracts are also interesting to consider regarding the proposed postero-anterior serial organization of frontal functions. A U-shaped organization of connections is concordant with caudal-rostral models [START_REF] Badre | Functional magnetic resonance imaging evidence for a hierarchical organization of the prefrontal cortex[END_REF][START_REF] Badre | Cognitive control, hierarchy, and the rostro-caudal organization of the frontal lobes[END_REF][START_REF] Christoff | Prefrontal organization of cognitive control according to levels of abstraction[END_REF]) of lateral prefrontal organization, such as the cascade model of cognitive control [START_REF] Koechlin | An information theoretical approach to prefrontal executive function[END_REF][START_REF] Koechlin | The architecture of cognitive control in the human prefrontal cortex[END_REF][START_REF] Koechlin | The role of the anterior prefrontal cortex in human cognition[END_REF]. In these models, more posterior regions correspond to a lower-level module and receive signals from a higherlevel module supported by more anterior regions. The processing-level of these modules may depend on the level of cognitive control, of relational complexity or of abstraction of representations. The most anterior region, or rostral prefrontal cortex, is thought to support the highest level modules, and has indeed been consistently shown involved in the highest level of cognitive abilities, such as relational reasoning [START_REF] Smith | Localizing the rostrolateral prefrontal cortex at the individual level[END_REF][START_REF] Vartanian | Dissociable neural systems for analogy and metaphor: implications for the neuroscience of creativity[END_REF][START_REF] Volle | Specialization of the rostral prefrontal cortex for distinct analogy processes[END_REF]), combination of remote information [START_REF] Gonen-Yaacovi | Rostral and caudal prefrontal contribution to creativity: a meta-analysis of functional imaging data[END_REF], coordinating goals and subgoals [START_REF] Koechlin | The role of the anterior prefrontal cortex in human cognition[END_REF][START_REF] Kroger | Recruitment of anterior dorsolateral prefrontal cortex in human reasoning: a parametric study of relational complexity[END_REF] and multitasking [START_REF] Burgess | Mesulam's frontal lobe mystery re-examined[END_REF][START_REF] Volle | The role of rostral prefrontal cortex in prospective memory: a voxel-based lesion study[END_REF]. The U-shaped tracts composing the frontal longitudinal system may constitute the anatomical link between successive modules. This hypothesis may be tested in future studies.

It is also interesting to note that projection areas of different tracts are not randomly distributed in the frontal lobes, but seem to overlap in preferred locations, such as the inferior frontal junction, the posterior part of the superior frontal sulcus, the inferior frontal gyrus, and the rostral prefrontal region. These cortical regions have been well observed in a myriad of functional imaging studies, using various cognitive tasks. They may well be functional hubs, in addition to anatomical hubs.

White matter correlations with age and education.

We report that age-related white matter differences distribute unevenly, inducing subtle differences in several tracts of the frontal region but not all. Our analysis confirms previous evidence of reduced microstructural organization in the frontal portion of the corpus callosum associated with aging [START_REF] Lebel | Longitudinal development of human brain wiring continues from childhood into adulthood[END_REF][START_REF] Lebel | Microstructural maturation of the human brain from childhood to adulthood[END_REF][START_REF] Lebel | Age-related regional variations of the corpus callosum identified by diffusion tensor tractography[END_REF]). This commissural decline may explain why increased reaction time is associated with aging, as reported during tasks requiring interhemispheric transfer [START_REF] Gootjes | Age effects in identifying and localising dichotic stimuli: a corpus callosum deficit?[END_REF][START_REF] Reuter-Lorenz | Differential effects of aging on the functions of the corpus callosum[END_REF] that may sometimes be compensated by alternative use of a different portion of the corpus callosum [START_REF] Schulte | Fiber tract-driven topographical mapping (FTTM) reveals microstructural relevance for interhemispheric visuomotor function in the aging brain[END_REF]. Similarly our analysis confirmed the decreased age-related microstructural measurements of the left and right corticospinal tract recently associated with age-related differences in fine motor control [START_REF] Holtrop | Investigating Age-related changes in fine motor control across different effectors and the impact of white matter integrity[END_REF].

Age-related differences reported in the anterior thalamic radiations bilaterally are concordant with the age-related significant decrease of the neuron number in the anterior thalamic nuclei [START_REF] Panadero | Memory and aging: changes in the mammillary body and anterior thalamic nuclei due to age][END_REF]. The anterior thalamic radiations are also part of network involving the hippocampus and dedicated to episodic memory [START_REF] Aggleton | Hippocampalanterior thalamic pathways for memory: uncovering a network of direct and indirect actions[END_REF][START_REF] Catani | A revised limbic system model for memory, emotion and behaviour[END_REF]. Hence the decreased age-related microstructural integrity of the anterior thalamic radiation might be related to the normal aging of episodic memory function [START_REF] Brickman | Aging and memory in humans[END_REF][START_REF] Charlton | The relationship between episodic long-term memory and white matter integrity in normal aging[END_REF].

We also explored age-related correlations for new pathways not defined in other atlases.

For the first time to our knowledge, we reported an age-related microstructural decreased integrity of two interlobar tracts in the right hemisphere (SLF I and III) and three intralobar tracts in the left hemisphere (the frontal inferior longitudinal fasciculus, the frontal orbito-polar tract and the fronto-insular tract 5). The contrast between the two hemispheres suggest that white matter differences related to age are not similar for the left and the right hemispheres and may be supported by different mechanisms.

Education level did not correlate significantly with tract-specific measurements, even when controlling for the effect of age. This finding suggests that education does not change or strengthen the white matter microstructure of the frontal lobe. Compensation mechanisms to counteract neural decline may therefore be supported by other mechanisms such as the recruitment of supplementary "compensatory" brain networks [START_REF] Stern | Brain networks associated with cognitive reserve in healthy young and old adults[END_REF][START_REF] Stern | A common neural network for cognitive reserve in verbal and object working memory in young but not old[END_REF][START_REF] Cabeza | Aging gracefully: compensatory brain activity in high-performing older adults[END_REF][START_REF] Cabeza | Frontal lobes and aging: Deterioration and compensation[END_REF]. Comparing the effect of education or measures of general intelligence using task-specific functional connectivity in young adults and elderly may provide further support to this hypothesis.

FA and HMOA tract specific measurements.

Spherical deconvolution tractography produced larger bundles, including voxels where fibers are crossing, which indubitably decreased FA and may have biased FA measures. New tract-specific measurements, such as HMOA, are less susceptible to crossing fiber effects.

Correlations between age and tract-specific HMOA, but not with FA, suggest that HMOA produces a better estimate of white matter bundle microstructure, when crossing with other fiber populations (Dell'Acqua et al. 2013). Therefore, when undertaking tract-specific measurements using advanced fiber-crossing algorithms, it is preferable to use real tract-specific measures, such as HMOA, rather than voxel specific measures, such as FA.

Limitations

We employed a structure-specific approach to quantify age-related changes occurring within frontal tracts instead of a whole brain voxel based approach. This approach may have missed differences that may occur at the voxel level. However structure-specific approach does not suffer normalization approximations and is therefore more sensitive than analysis performed over the whole brain [START_REF] Catani | Diffusion tensor magnetic resonance imaging tractography in cognitive disorders[END_REF]). New hybrids approaches, such as the one reported in Yushkevich et al. (Yushkevich et al. 2008), using true-tract specific indices to characterize white matter diffusion when crossing might bring the best of structure specific and voxel based approaches.

We noted that volume measurements have a larger standard deviation than FA or HMOA measurements. This is important as a larger standard deviation reduces the statistical power and decreases the chances of reporting a significant result. Accuracy of tract volume measurements depends on the quality of the dissection, which relies on the reliability of the operator and the quality of the dataset. We employed a semi-automatic approach for the dissection in order to reduce variability, and subsequently, all tracts were visually inspected to reduce false positives.

However false negatives may still have been the source for this increase in volume variability [START_REF] Jones | White matter integrity, fiber count, and other fallacies: the do's and don'ts of diffusion MRI[END_REF]Kristo et al. 2013a).

We also report preliminary results on the distribution of the corpus callosum, corticospinal tract, striatum, pons and thalamus connections to the cortical surface. Axonal tracing studies suggest that a comparable distribution of projection for the corpus callosum [START_REF] Myers | Organization of forebrain commissures[END_REF]), the striatum [START_REF] Haber | Striatonigrostriatal pathways in primates form an ascending spiral from the shell to the dorsolateral striatum[END_REF], the thalamus (Barbas et al., 1991) and the pons [START_REF] Schmahmann | Prefrontal cortex projections to the basilar pons in rhesus monkey: implications for the cerebellar contribution to higher function[END_REF] also exists in monkeys. Indeed, the frontal corpus callosum decreases its probability of projection for the frontal eye-field, the orbito frontal gyrus and the frontal operculum in monkey axonal tracing [START_REF] Myers | Organization of forebrain commissures[END_REF] as well as in the map we presented.

Similarly, the striatal projection seems to be reduced for the Brodmann areas 9 and 46 [START_REF] Haber | Striatonigrostriatal pathways in primates form an ascending spiral from the shell to the dorsolateral striatum[END_REF] and lighter projections from the pons were received from the medial and ventrolateral cortices [START_REF] Schmahmann | Prefrontal cortex projections to the basilar pons in rhesus monkey: implications for the cerebellar contribution to higher function[END_REF]. However, quantitative maps of the fiber density arising from or terminating in the simian cortex are missing in literature and prevent us from drawing firm comparative anatomy conclusions. Furthermore, in most of our participants the corticospinal tract did not connect to the whole motor homunculus. Hence, although spherical deconvolution may put forward an encouraging improvement in the quality of the dissections as compared to standard tensor tractography, it is not exempt of limitations. It should be used with caution, especially when reconstructing projection pathways or in clinical applications where edema, tumors, lesions and infarcts are present then multiple compartments and partial volume effects will limit and confound spherical deconvolution [START_REF] Roine | Isotropic non-white matter partial volume effects in constrained spherical deconvolution[END_REF].

Finally, our cross sectional report suggests heterogeneous age-related microstructural differences within the frontal lobe, that ideally, will require replication in a longitudinal study.

Conclusions

This atlas of the frontal lobe connections may help identifying tracts connecting functional activations thus contributing to a better understanding of frontal functions and processes. It may also help clinicians to diagnose the disconnection or abnormality of specific tracts on patients' MRIs following a wide range of pathologies including but not limited to stroke, degenerative diseases, traumatic brain injury, multiple sclerosis, tumors, vascular malformations, leukoencephalopathies, and infectious diseases. In relation to clinical disorders, our maps could increase understanding of cognitive deficits, as well as mechanisms of brain recovery and plasticity. However, special caution should be taken when overlapping the atlas maps with the structural data of patients due, for example, to tissue displacement and an altered anatomy [START_REF] Brett | Spatial normalization of brain images with focal lesions using cost function masking[END_REF][START_REF] Clark | White matter fiber tracking in patients with space-occupying lesions of the brain: a new technique for neurosurgical planning?[END_REF]. Finally, tract-specific correlations with aging produce interesting hypotheses on possible early cognitive decline dissociations that may be tested in stroke and elderly patients.

West RL (1996)

Tract specific measurements

de

  Schotten et al. 2014a) or disconnections caused by strokes (Thiebaut deSchotten et al. 2014b).

Fig. 2 .

 2 Fig. 2. Cortical projection in the MNI stereotaxic space of the frontal interlobar tracts [top to

Fig. 3 .

 3 Fig. 3. Cortical projection in the MNI stereotaxic space of the commissural and projection tracts.

Fig. 4 .

 4 Fig. 4. Cortical projection in the MNI stereotaxic space of frontal short U-shaped tracts (top to

Fig. 5 .

 5 Fig. 5. Correlation between the participants' age and tract-specific measurements. Participants

Fig. 6 .

 6 Fig. 6. Correlation between the participants' age and tract-specific measurements. Participants

  An application of prefrontal cortex function theory to cognitive aging. Psychological bulletin 120 (2):272-292. Xiong YY, Mok V (2011) Age-related white matter changes. J Aging Res 2011:617927. Yushkevich PA, Zhang H, Simon TJ, Gee JC (2008) Structure-specific statistical mapping of white matter tracts. Neuroimage 41 (2):448-461.

	CAPTIONS
	Fig. 1. Spherical deconvolution tractography of frontal white matter connections. For each tract a
	single participant map is supplied as a representative example of the individual anatomy. First
	(SLF I, light blue, n°1), second (SLF II, dark blue, n° 2) and third (SLF III, pink, n°3) branches
	of the superior longitudinal fasciculus; cingulum (yellow, n°4); uncinate (pink, n°5); long (LS,
	red, n°6) and anterior (AS, green, n°7) segments of the arcuate fasciculus; inferior fronto-
	occipital fasciculus (IFOF, n°8); frontal corpus callosum (CC, n°9); corticospinal tract (CST,
	n°10); fronto-thalamic projections or anterior thalamic radiations (blue, n°11); fronto-striatal
	projections (yellow, n°12); fronto-pontine projections (n°13); paracentral U tract (pink, n°14);
	hand superior (green, n°15), middle (yellow, n°16) and inferior (red, n°17) U tract; face U tract
	(blue, n°18); frontal aslant tract (pink, n°19); fronto insular tract 1 (FIT 1, light blue, n°20), 2
	(FIT 2, dark blue, n°21), 3 (FIT 3, red, n°22), 4 (FIT 4, yellow, n°23) and 5 (FIT 5, green, n°24);
	frontal superior longitudinal (FSL, light blue n°25); frontal inferior longitudinal (FIL, pink,
	n°26), frontal orbito-polar tract (FOP, yellow, n°27), fronto-marginal tract (FMT, red, n°28).

ACKNOWLEDGMENTS

We thank Dr. Marco Catani for his assistance with the anatomical dissection and insightful discussions. We also thank the French Agence Nationale de la Recherche for their support of this project (project CAFORPFC, no. ANR-09-RPDOC-004-01 and project PHENOTYPES, no.

ANR-13-JSV4-0001-01). In addition, we also thank the program "Investissements d'avenir" (ANR-10-IAIHU-06) for their generous support.

SUPPLEMENTARY MATERIAL

Atlasing the frontal lobe connections and their variability due to age and education: a spherical deconvolution tractography study.

Rojkova K. 1,2 , Volle E. 1 , Urbanski M.

DIFFUSION-WEIGHTED IMAGING ACQUISITION AND PREPROCESSING

A total of 70 near--axial slices were acquired on a Siemens 3 Tesla VERIO TIM system equipped with a 32--channel head coil. We used a fully optimised acquisition sequence for the tractography of diffusion--weighted imaging (DWI), which provided isotropic (2 × 2 × 2 mm) resolution and coverage of the whole head with a posterior--anterior phase of acquisition. The acquisition was peripherally--gated to the cardiac cycle [START_REF] Conturo | Diffusion MRI: precision, accuracy and flow effects[END_REF][START_REF] Jones | Spatial normalization and averaging of diffusion tensor MRI data sets[END_REF][START_REF] Turner | Echo--planar imaging of intravoxel incoherent motion[END_REF]) with an echo time (TE) = 85 msec. We used a repetition time (TR) equivalent to 24 RR. At each slice location, 6 images were acquired with no diffusion gradient applied. Additionally, 60 diffusion--weighted images were acquired, in which gradient directions were uniformly distributed on the hemisphere with electrostatic repulsion. The diffusion weighting was equal to a b--value of 1500 sec mm -2 . One supplementary image with no diffusion gradient applied but with reversed phase--encode blips was collected. This provides us with a pair of images with no diffusion gradient applied with distortions going in opposite directions. From these pairs the susceptibility--induced off--resonance field was estimated using a method similar to that described in (Andersson et al. 2003) and corrected on the whole diffusion weighted dataset using the tool TOPUP as implemented in FSL [START_REF] Smith | Advances in functional and structural MR image analysis and implementation as FSL[END_REF]. Finally, at each slice, diffusion-weighted data were simultaneously registered and corrected for subject motion and geometrical distortion adjusting the gradient accordingly (ExploreDTI http://www.exploredti.com) see [START_REF] Leemans | The B--matrix must be rotated when correcting for subject motion in DTI data[END_REF]

SPHERICAL DECONVOLUTION TRACTOGRAPHY RECONSTRUCTION

Spherical deconvolution was chosen to estimate multiple orientations in voxels containing different populations of crossing fibers (Alexander 2006;[START_REF] Stern | Brain networks associated with cognitive reserve in healthy young and old adults[END_REF]Tournier et al. 2004). The damped version of the Richardson--Lucy algorithm for spherical deconvolution (Dell'acqua et al. 2010) was calculated using a software that was developed in--house an free on demand to the authors. Algorithm parameters were chosen as previously described (Dell'Acqua et al. 2013). A fixed--fiber response corresponding to a shape factor of α = 1.5 × 10-3 mm2/s was chosen (Dell'Acqua et al. 2013). For each participant, a convergence speed (CS) map of the deconvolution algorithm [START_REF] Dell'acqua | Convergence Maps from Richardson--Lucy Spherical Deconvolution Algorithm for the Detection of White Matter in HARDI[END_REF]) was estimated. CS map quantifies how quickly the residual fitting error between the diffusion signal, and the fiber model as identified by the deconvolution algorithm decays within each voxel. CS maps better contrast white matter regions showing a smaller partial volume effect, as compared to FA or similar anisotropy maps. Fiber orientation estimates were obtained by selecting the orientation corresponding to the peaks (local maxima) of the fiber orientation distribution (FOD) profiles. To exclude spurious local maxima, we applied both an absolute and a relative threshold on the FOD amplitude. A first "absolute" threshold was used to exclude intrinsically small local maxima due to noise or isotropic tissue. This threshold was set to 3 times the mean amplitude of a spherical FOD obtained from a gray matter isotropic voxel (and therefore also higher than an isotropic voxel in the cerebro--spinal fluid). A second "relative" threshold of 10% of the maximum amplitude of the FOD was applied to remove the remaining local maxima with values greater than the absolute threshold [START_REF] Dell'acqua | Mapping Crossing Fibres of the Human Brain with Spherical Deconvolution: Towards an Atlas for Clinico--Anatomical Correlation Studies Proceedings of the International Society of Magnetic Resonnance[END_REF]).

Whole brain tractography was performed selecting every brain voxel with at least one fiber orientation as a seed voxel. From these voxels, and for each fiber orientation, streamlines were propagated using Euler integration with a step size of 1 mm (as described in Dell' Acqua et al. 2013). When entering a region with crossing white matter bundles, the algorithm followed the orientation vector of least curvature (as described in Schmahmann and Pandya 2007) Streamlines were halted when a voxel without fiber orientation was reached or when the curvature between two steps exceeded a threshold of 60°. Spherical deconvolution, fiber orientation vector estimations and tractography were performed using in--house software developed with Matlab 7.8 (http://www.mathworks.com).

TRACTOGRAPHY DISSECTIONS

We produced a CS--MNI152 template using Advanced Normalization Tools (ANTs, http://stnava.github.io/ANTs/) (Avants and Gee 2004;[START_REF] Avants | The optimal template effect in hippocampus studies of diseased populations[END_REF], which build a template iteratively combining affine and diffeomorphic deformations (Avants et al. 2008;[START_REF] Klein | Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration[END_REF]. Participants were first realigned to the MNI using affine transformation. Then the template was built from all the subjects iteratively (n=4) using diffeomorphic deformations.

In order to facilitate the tractography dissection, regions of interest (ROI) based on the guidelines provided in previous reports (anterior thalamic projections see [START_REF] Behrens | Non--invasive mapping of connections between human thalamus and cortex using diffusion imaging[END_REF], main association fibers see (Catani and Thiebaut de Schotten 2008); subdivision of the arcuate see (Catani et al. 2007); fronto--pontine fibers see (Catani et al. 2008); frontal U--shape fibers see (Catani et al. 2012); fronto--striatal projections see [START_REF] Cohen | Connectivity--based segregation of the human striatum predicts personality characteristics[END_REF]; SLF subdivision see (Thiebaut de Schotten et al. 2011a)) were defined on the CS--MNI template calculated above. In addition we performed tractography reconstruction of all the streamlines emerging from the frontal lobe and all streamlines not involving the frontal lobes. Detailed guidelines about the dissection can be found in supplementary online methods. For each participant, the CS Map was registered to the CS--MNI152 template using ANTs. The inverse deformation was then applied to the ROIs, as defined in the CS--MNI152, in order to bring them to the native space of every participant Individual dissections of the tracts were then visually inspected and corrected in trackvis (http://www.trackvis.org), when required, by two anatomists (KR and MTS) who manually readjusted the ROIs using rigid transformation. Tract--specific measurements, such as the relative volume (volume of the tract divided by the total volume of the brain), the fractional anisotropy (mean FA value for the whole tract, [START_REF] Basser | MR diffusion tensor spectroscopy and imaging[END_REF] and a new index-- employed as a surrogate for tract microstructural organization (i.e. mean Hindrance Modulated Orientational Anisotropy or HMOA for the whole tract Dell'Acqua et al. 2013) were extracted from each dissected tract. HMOA has the advantage of providing information about the microstructural diffusion properties of distinct fiber orientations and therefore specific to the orientation of the reconstructed tracts and more accurate than classical fractional anisotropy measures, which decreases when fibers cross due to local partial volume effect. 1 -Frontal and non-frontal tracts Supplementary Fig. 1: Delineation of the ROI used for the tractography of the frontal tracts and the non--frontal tracts. Left bullet, color legend for the ROI; right bullet, Boolean rules employed for the dissection. A one--ROI approach was used to isolate the frontal lobe streamlines. We draw a region of interest including the whole frontal lobe cortical ribbon. We only displayed streamlines that were either emerging or ending inside the frontal lobe cortical ribbon. For the non--frontal tracts we excluded from the whole brain tractography the streamlines connecting the frontal lobe cortical ribbon. A multiple region of interests (ROIs) approach was used to isolate the three components of the superior longitudinal fasciculus (i.e. SLF I, SLF II and SLF III). Three coronal ROIs were delineated around the white matter of the superior, middle and inferior/precentral frontal gyri, and another three 'AND' ROIs were delineated posteriorly in the parietal region. Streamlines of the arcuate fasciculus projecting to the temporal lobe were excluded using an axial 'NOT' ROI in the temporal white matter (the arcuate is not part of the longitudinal system as it forms an arc to reach the temporal lobe) 3 -Cingulum Supplementary Fig. 3 A two ROIs approach was used to isolate the frontal projections of the corpus callosum and the left and right fronto--striatal, anterior thalamic radiation and fronto--pontine fibers. A 'AND' sagittal ROIs was delineated around the white matter of the corpus callosum. We also defined three 'AND' region around the whole striatum (including the caudate and the putamen), thalamus, pons whole and frontal lobe cortical ribbon and only displayed streamlines that were either emerging or ending inside these ROIs. A two ROIs approach was used to isolate the paracentral, hand superior, hand middle, hand inferior and face U--shaped in the left and the right hemispheres. Eight 'AND' axial ROIs were delineated in both hemispheres around the white matter of the anterior and posterior portion of the paracentral gyrus, the superior portion of the omega sign designating the motor hand area and its corresponding post--central area, the inferior portion of the omega sign and its corresponding post--central area, the most inferior portion of the pre--and post-central gyri. Note that each axial ROI was delineated on three consecutive slices (approximately 6 millimeters thick). An upper length threshold of 50 millimeters was applied to the dissection to discard spurious reconstructions. 9 -Frontal aslant tract, superior and inferior frontal longitudinal fasciculi Supplementary Fig. 9: Delineation of the ROIs used for the tractography of the left and right frontal aslant tract, superior and inferior frontal longitudinal fasciculi. Left bullet, color legend for the 'AND' ROIs; right bullet, Boolean rules employed for the dissection. A two ROIs approach was employed to isolate the left and right frontal aslant tract, superior and inferior frontal longitudinal fasciculi. For each hemisphere we delineated an axial 'AND' ROI around the white matter of the superior frontal gyrus and a sagittal 'AND' ROI around the white matter of the inferior frontal gyrus (also including the pars opercularis, triangularis and orbitalis). Finally a 'AND' ROI was delineated around the white matter of the middle frontal gyrus in axial and sagittal sections. 10 -Fronto-insular tracts Supplementary Fig. 10: Delineation of the ROIs used for the tractography of the left and right fronto--insular tracts. Left and middle bullets, color legend for the 'AND' ROIs; right bullet, Boolean rules employed for the dissection. A two ROIs approach was used to isolate the fronto--insular tracts. Five sagittal 'AND' ROI were delineated in both hemisphere around the white matter of the pars orbitalis, pars triangularis, pars opercularis, inferior portion of the precentral gyrus, and the subcentral gyrus. We also defined another 'AND' ROI around the whole insula displayed streamlines that were either emerging or ending inside this ROI. An upper length threshold of 50 millimeters was applied to the dissection to discard spurious reconstructions. Supplementary Fig. 11: Delineation of the ROIs used for the tractography of the left and righ frontomarginal and orbitopolar and tracts. Left bullet, color legend for the 'AND' ROIs; right bullet, Boolean rules employed for the dissection. A two--ROIs approach was used to isolate the frontal orbitopolar tract. Two coronal 'AND' ROI were delineated in both frontal lobes around the white matter of the lateral and the medial frontal pole separated by the frontomarginal sulcus, pars triangularis, pars opercularis, inferior portion of the precentral gyrus. For the orbitopolar tract we also defined another 'AND' ROI around the white matter of the posterior orbital gyrus.

-Superior longitudinal fasciculus Supplementary

MAPPING OF THE FRONTAL LOBE CONNECTIONS

Binary individual visitation maps were created for each tract by assigning each voxel a value of 1 or 0, depending on whether the voxel was intersected by the streamlines of the tract. Binary visitation maps of each dissected tracts were normalized to the MNI space using the same affine and diffeomorphic deformations as calculated above. We created percentage overlap maps by adding the normalized visitation maps from each subject at each point in the MNI space. Therefore, the overlap of the visitation maps varies according to inter-subject variability. A similar approach has been reported in Thiebaut de Schotten et al. (Thiebaut de Schotten et al. 2011b) 3D brain rendering of the MNI152 template provided with the FSL package was calculated using MRIcroGL (http://www.cabiatl.com/mricrogl/). Percentage overlap maps of the tracts were projected onto the 3D rendering of the MNI152, in order to visualize cortical projections (Thiebaut de Schotten et al. 2011a). Probabilistic maps from the Jubrain (https://www.jubrain.fz--juelich.de) were employed to identify the corresponding Brodmann areas. The overall visualization and screenshots were performed in MRIcroGL.

AGE-RELATED SHAPE CHANGES

In this section we assessed whether the shape of the tracts change over age. We assessed this question by producing tractography maps for each decade (Supplementary Fig. 12b). We also quantified similarities using cross-correlation between the maps for the 21--30 year--old maps and the maps for the other decades (Supplementary Fig. 12a). The result indicates that tracts reconstructed with tractography show a very high anatomical correspondence between decades and no decrease of this correspondence with age. Hence tractography in our sample did no argue for age related changes in the shape of the tract. Supplementary Fig. 12