Frédéric Desprez

George S Markomanolis

Frédéric Suter

Improving the Accuracy and Efficiency of Time-Independent Trace Replay

Keywords: Performance prediction, MPI, Simulation Prédiction de performances, MPI, Simulation

Simulation is a popular approach to obtain objective performance indicators on platforms that are not at one's disposal. It may help the dimensioning of compute clusters in large computing centers. In a previous work, we proposed a framework for the off-line simulation of MPI applications. Its main originality with regard to the literature is to rely on time-independent execution traces. This allows us to completely decouple the acquisition process from the actual replay of the traces in a simulation context. Then we are able to acquire traces for large application instances without being limited to an execution on a single compute cluster. Finally our framework is built on top of a scalable, fast, and validated simulation kernel. In this paper, we detail the performance issues that we encountered with the first implementation of our trace replay framework. We propose several modifications to address these issues and analyze their impact. Results shows a clear improvement on the accuracy and efficiency with regard to the initial implementation.

Introduction

Understanding the performance and behavior of parallel applications relying on the Message Passing Interface (MPI) on various compute infrastructures is an main concern for many developers. Many tools exist, such as TAU [START_REF] Shende | The Tau Parallel Performance System[END_REF], Scalasca [START_REF] Geimer | A Scalable Tool Architecture for Diagnosing Wait States in Massively Parallel Applications[END_REF], or the joint effort on SCORE-P [START_REF] Mey | Score-P -A Unified Performance Measurement System for Petascale Applications[END_REF], that provide strong analysis methods based on profiling to detect performance bottlenecks.

Simulation is another popular approach for predicting the performance of an application. When a platform is yet to be specified and purchased, simulations can be used to determine a cost-effective hardware configuration appropriate for the expected application workload. Conversely, simulations can also be used to study the performance behavior of an application by varying the hardware characteristics of an hypothetical platform. In a classroom setting, students without access to a parallel platform could execute applications in simulation on a single node as a way to learn the principles of parallel programming and high-performance computing. Simulation of an application on a platform may also be useful even when the platform is available. For instance, the simulation may bypass actual computations performed by the application, and only simulate the corresponding delays of these computations. In this case the simulated application produces erroneous results, but its performance behavior may be preserved. It is then possible to conduct development activities for performance tuning in simulation only on a small-scale platform, thereby saving time when compared to real executions on a large-scale platform. Furthermore, access to large-scale platforms is typically costly (possible access charges to the user, electrical power consumption). The use of simulation can thus not only save time but also money and resources.

Many simulation frameworks have been proposed over the last decade. They fall into two categories: on-line simulation, also called simulation via direct execution, and off-line simulation, also called postmortem simulation. In on-line simulation the application is executed but part of the execution takes place within a simulation component. In off-line simulation a trace of a previous execution of the application is "replayed" on a simulated platform. Most of the existing tools in this category rely on timed traces, i.e., each traced event is associated to a time-stamp. This approach creates a tight link between the trace and the set of machines used to produce it. This link is a clear limit to the applicability of the offline approach. In [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF] we proposed a first prototype of a framework that breaks this link by getting rid of the time-stamps. This greater flexibility related to the acquisition of traces comes with some issues of its own. Some were highlighted by the experiments conducted in [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF]. First, the acquisition of the traces is sometimes time-consuming and may, in very specific cases, be close to the execution time of the application itself. Then, as our framework relies solely on volumes and not on timestamps, it requires a very accurate measure of hardware counter values, and a minimal impact of the instrumentation. We noted a certain discrepancy of the measured number of instructions in our first implementation that may hinder the quality of the simulation results. Finally, while the accuracy of the simulations is very good for either compute or communication intensive applications, it was not the case for applications that alternate computations and communications.

This paper is organized as follows. In Section 2 we detail the issues raised by the first prototype of our off-line simulation framework based on Time-Independent Traces. In Section 3 we present the different modifications of the framework aiming at addressing these issues. The impact of the proposed modifications are given in Section 4. We review related work in Section 5. Finally, we draw some conclusions and detail future work directions in Section 6.

Identifying the Issues

In this section, we analyze the performance of our Time-Independent Trace replay framework. We detail the different issues that raised in the evaluation section of [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF]. We reproduce experiments conducted on the same compute cluster, called bordereau, that comprises 93 2.6 GHz Dual-Proc, Dual-Core AMD Opteron 2218 nodes interconnected through a single 10 Gigabit switch. Each core has a L2 cache of 1 MB. As this cluster is now aging and is prone to failures and suspect behaviors, we complement this study with results obtained on a more recent cluster, called graphene, that comprises 144 2.53GHz Quad-Core Intel Xeon x3440 nodes. The L2 cache owned by each core is two times larger than on the bordereau cluster. The nodes are scattered across four cabinets, and interconnected by a hierarchy of 10 Gigabit switches.

Instrumentation Time Overhead

In our first implementation of a Time-Independent Trace replay framework, we relied on the TAU profiling tool [START_REF] Shende | The Tau Parallel Performance System[END_REF] and its Program Database toolkit (PDT) [START_REF] Lindlan | A Tool Framework for Static and Dynamic Analysis of Object-Oriented Software with Templates[END_REF] to trace the execution of an MPI application. We also took advantage of the selective instrumentation feature offered by TAU to focus on specific parts of the studied application. Using such tools was very convenient as the instrumentation is done automatically. However, the studied application is instrumented at a very fine grain. This may lead to a non-negligible, overhead on the execution time as shown in [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF]. Basically, this overhead is the difference between the respective execution times of an application with and without instrumentation. As identified in [START_REF] Mussler | Reducing the overhead of direct application instrumentation using prior static analysis[END_REF], this overhead has different sources, e.g., intrusion of the measurement mechanisms, building the complete call path of the application and flushing the trace on disk.

In [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF], we measured an overhead ranging from 4.1% to 39.5%, depending on the problem instance, for the execution of a LU factorization. The largest value was observed for small problem sizes executed on a large number of compute resources. While this overhead remains acceptable as traces are acquired once to be replayed multiple times, we aim at reducing it as much as possible. On the graphene cluster, this overhead ranges from 9.8% to 25.5% in similar conditions. For 128 processes, the overhead increases to 37.9%.

Impact of Instrumentation on Hardware Counters Values

A second issue is related to the impact of the instrumentation on the values of the hardware performance counters used to produce a Time-Independent Trace. In our work, we use the total number of instructions executed by the studied application. Adding instrumentation probes may increase this value and then impact the accuracy of a simulation based on the produced trace. To measure the impact of a fine grain instrumentation such as that done by TAU, we perform the following simple experiment. We execute two versions of the same application, namely the LU factorization, on the same compute cluster. The first version is fully instrumented using TAU. In the second version, we just insert calls to get the value of the hardware performance counter that measures the number of executed instructions at the beginning and end of the studied section of the code. Figure 1 shows the distribution across processes of the relative difference in terms of measured numbers of executed instructions between fine and coarse grain instrumented versions of various instances of the LU factorization. The comparison between the two executions is done on a per-process basis. Moreover, we ran ten runs of each version and display the average values. Finally the average total number of instructions per process measured with the coarse grain instrumentation method is 1.70e11 for the smallest instance (B-8) while it is 8.87e10 for the largest instance (C-64).These results were obtained on the bordereau cluster.

With the exception of instance B-64, we see that the fine grain instrumentation leads to an increase of the measured number of instructions of around 10-13%. As our Time-Independent Traces are directly using these numbers of instructions, such an increase will have a negative impact on the accuracy of the simulation. In other words, it will likely simulate something closer to the instrumented version than the original application. The discrepancy becomes even worse when communications start to prevail. This is the case with the B-64 instance in which the balance between computation and computation changes due to data distribution. Figure 2 presents similar results but obtained on the more recent graphene cluster. For instances that involves between 8 and 64 processes, we note a similar discrepancy on this cluster as on the bordereau cluster. It ranges from 11 to 16%. However, we observe a clear increase of the discrepancy as the number of processes grows for the Class B instances. For 128 processes, it raises up to 23%. It can easily be explained as, for such an instance, the rather small input data is highly distributed. Each process thus only has a really small amount of data to compute on and is more impacted by the instrumentation.

Calibration and Cache Usage

An essential step to make accurate performance predictions through trace replay is the calibration of the simulation framework. In our framework, it consists in determining the number of instructions a CPU can compute in one second and the latency and bandwidth of communication links. Such a calibration strongly depends on both application and execution environment.

In our first implementation, we chose to do the calibration based on the execution of a small instance of the studied application. For instance, we consider the A-4 instance to calibrate our tool for the simulation of larger instances of the LU factorization. The problem there is that the A-4 instance is small enough to have all data stored in the L2 cache. Then it does not allow us to capture the worse performance achieved by instances for which the data distribution leads to data movements between L2 cache and main memory.

The instruction rate being one the most important parameter to achieve accurate simulated executions, we have to define a more complex calibration step to better take cache usage in account.

Accuracy of the Simulated Execution

In the evaluation of the performance of our off-line simulation framework we conducted in [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF], we compared the execution time of various instances of the LU factorization to the simulated time obtained by replaying Time-Independent Traces (see Fig. 8 in [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF]). For such an application that is composed of millions of interleaved computation and communication actions, the accuracy of the simulated replay was not as good as expected. Figure 3 shows the evolution of the relative error between the time to compute a LU factorization and its simulated counterpart when the number of involved processes vary. We see that the inaccuracy of the simulated version is not stable, but increases rather linearly with the number of processes. Moreover, our tool underestimates the execution time for small number of processes, i.e., when each process owns a larger share of data, and then overestimates the execution as the number of processes grows.

This inaccuracy may have different origins related either to computations or communications. Such an inaccuracy may for instance come from the discrepancy of the measured number of instructions that we just mentioned. Indeed it directly impacts the calibration of the replay tool that determines the rate at which each machine can process instructions. Among other potential sources, we assume that every part of the application can be processed at the same rate which may not be the case for some applications, due to cache affinity for instance. However, the general trend shown by Figure 3 indicates that the main source of inaccuracy comes from a bad estimation of communication times. The most probable source might be that the LU factorization, as implemented in the NAS Parallel Benchmark suite, implies a lot of small messages (less than 64 KiB). Then most of the point-to-point communications use the eager mode of MPI in which the sender copies the data directly into the memory of the receiver without having to wait for it to post the corresponding receive. This particular communication protocol was implemented in a too simple way, based on asynchronous communications for small messages, in the first version of our framework. As the number of small messages increases with the number of involved processes, small inaccuracies tend to accumulate and lead to a large overall relative error.

Addressing the Issues

In the previous section we have highlighted four issues that all affect the performance and the quality of our off-line simulation framework. We detail hereafter four approaches to address these issues. The first lever we use is external to our framework, i.e., leveraging compiler optimizations. The other three imply modifications of the first implementation presented in [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF], i.e., reducing the intrusiveness of the instrumentation method, changing the simulation back-end to better handle communications, and improving the calibration procedure.

Compiler Optimizations

The first modification we made to our trace acquisition procedure is to activate compiler optimizations, typically by using the -O3 flag. Turning on an optimization flag makes the compiler attempt to improve the performance and/or code size at the expense of compilation time. Among the optimizations that may help to reduce the discrepancy in the measured number of instructions are the loop unrolling, vectorization, and function inlining.

Less Intrusive Instrumentation Method

As stated in Section 2.1, the default instrumentation provided by TAU leads to some execution time overhead. The main source of this overhead comes from the building of the complete call path of the instrumented application. While this information may be useful to users aiming at identifying performance bottlenecks in their application, it is useless to our specific usage of the produced traces. Indeed, the only information required by our Time-Independent Trace replay framework are: (i) the volume of computation, in number of instructions executed between two MPI calls; (ii) the name and parameters of each MPI call; and (ii) the volume of data transferred by each communication operation. This is illustrated by the trace snippet below. The selective instrumentation feature provided by TAU offers us a way to obtain all this required information, and only it. Then we modify our instrumentation method as follows. First we create a file, named exclude.pdt, in which we indicate to TAU the list of source files that have to be excluded from the detailed instrumentation. Actually, we exclude all the source files using the special character * .

BEGIN_FILE_EXCLUDE_LIST * END_FILE_EXCLUDE_LIST

Then, we add the following option to the command line to take this exclusion file into account.

-optTauSelectFile=/path/exclude.pdt Thanks to this straightforward modification, the instrumentation becomes minimal with regard to our specific needs. Indeed the performance hardware counter that measures the number of executed instructions will be triggered when entering and exiting MPI functions. All the information related to MPI calls, i.e., id of the process that calls this function and the function name and calling parameters, are traced exactly as in the previous implementation.

We expect this simple modification to have an impact both on the execution time overhead and the discrepancy of measured numbers of instructions between non-instrumented and instrumented versions of an application.

Simulation Back-end

Our simulation framework is is tightly connected to the SIMGrid project [START_REF] Casanova | SimGrid: a Generic Framework for Large-Scale Distributed Experiments[END_REF] that provides core functionalities for the simulation of distributed applications in heterogeneous distributed environments. SIMGrid relies on a scalable and extensible simulation engine and offers several user APIs. Our first implementation of a trace replay mechanism was based on the MSG API that was designed to implement simulators using Concurrent Sequential Processes (CSP). This choice forced us to mimic the behavior of MPI primitives, e.g., collectives operations or protocols depending on the message size, with some crude simplifications. More importantly it decouples this effort on the off-line simulation of MPI applications to that on the on-line simulation of MPI applications also taking place within the SIMGrid project.

In [START_REF] Clauss | Single Node On-Line Simulation of MPI Applications with SMPI[END_REF], we proposed a new simulator, SMPI, for the on-line approach. It is, as the present off-line approach, built on top of the simulation kernel of the SIMGrid toolkit and thus also benefits of its fast, scalable, and validated network models. SMPI extends the existing models with an original piece-wise linear model to take into account the specifics of the cluster interconnect. A salient property of SMPI is its capacity to simulate MPI applications on a single node.

We have reimplemented the trace replay mechanism directly within SMPI, instead of on top of the MSG API. This rewriting leads to many beneficial changes. First, it simplifies the way actions, i.e., events stored in the trace are passed to the simulation kernel. To illustrate this, we show the evolution of the implementation of the action that corresponds to an MPI_Send call. The format of this action is <id> send <dst_id> <volume>. static void action_send(const char * const * action){ char to[250]; double size=parse_double(action [START_REF] Casanova | SimGrid: a Generic Framework for Large-Scale Distributed Experiments[END_REF]); sprintf(to, "%s_%s", MSG_process_get_name(MSG_process_self()), action [START_REF] Badia | Dimemas: Predicting MPI applications behavior in Grid environments[END_REF]); if (size<65536) action_Isend(action); else MSG_task_send(MSG_task_create(NULL, 0, size, NULL), to); }

Inria

The underlying principle there is to create a mailbox named after the ids of the sender and receiver processes, create a task that corresponds to the communication of a message, and then send this task into that mailbox. A matching action on the receiver side, will read the contents of the mailbox and execute the task, which actually starts the simulated communication. This implementation does not reflect some the MPI specifics, especially for small messages. Indeed, when the message is smaller than 64KB, the eager mode is activated. This means that the sender does not have to wait for the receiver to allocate memory and directly copies the data remotely. We tried to model that by using an asynchronous send for such small messages. However, it is not what is actually implemented by most MPI runtimes.

Through the analyze of execution traces, we found that an MPI_send for small messages acts in a sort of "detached" way. From the application point of view, the send corresponds to the time of a copy of the data in the memory. Moreover, if the receive is issued after the send, the data is already stored in memory. There again, the application only sees the duration of a memory copy. Modeling this particular behavior has been done in SMPI. Doing the same in our initial trace replay tool would have been redundant, error prone, and above all difficult with the MSG API. The new implementation allows us to benefit of all the complex optimizations offered by SMPI in a seamless way. Indeed the send action is now: static void action_send(const char * const * action){ int to = atoi(action [START_REF] Badia | Dimemas: Predicting MPI applications behavior in Grid environments[END_REF]); double size=parse_double(action [START_REF] Casanova | SimGrid: a Generic Framework for Large-Scale Distributed Experiments[END_REF]); smpi_mpi_send(NULL, size, MPI_BYTE, to, 0, MPI_COMM_WORLD); }

The second main change implied by this rewriting is the user view of the replay framework. In the former version everything was exposed to the user. Everything is now embedded and replaying a Time-Independent Trace with SimGrid simply amounts to run the following program.

int main(int argc, char * argv[]){ smpi_replay_init(&argc, &argv); smpi_action_trace_run(NULL); smpi_replay_finalize(); return 0; } This code, that initializes some data structures, loads a trace, and destroys the data structures, is considered as a regular SMPI program. Then it is compiled with the smpicc wrapper and launched thanks to the smpirun command (see [START_REF] Clauss | Single Node On-Line Simulation of MPI Applications with SMPI[END_REF] for details), as follow smpirun -np 8 -hostfile hostfile -platform platform.xml \ ./smpi_replay trace_description where np and hostfile are classical parameters of mpirun. In addition, SMPI requires a description of the simulated platform, given in the platform.xml file. Our replay tool needs a single parameter, a file that list the names of the trace files to associate to each process. If this file contains a single entry, all the processes will look for the actions they have to perform into the same trace.

The last important change is about the format of the recv action. To keep our implementation as simple as possible, we had to add the message size to the parameters of this action. As we were already able to track this information in the traces produced by TAU, it required only a small modification of the existing extraction script.

Cache-Aware Calibration Step

As mentioned in Section 2.3, the first implementation of our framework relied on small instances, typically Class A on 4 processes for the NAS Parallel Benchmark suite, to calibrate the instruction rate of the simulated platform. The rationale was to limit the number of resources and the time needed to do the calibration as much as possible. Experiments shown that while using only as few resources as four cores did not raise any issue, limiting the calibration to a small problem size did hide some important phenomena.

In the particular case of the LU factorization, as soon as the share of the matrix owned by each process exceeds the capacity of the L2 cache, the performance drops, with a direct impact on the instruction rate. To circumvent this issue, we propose to complete our calibration procedure, when it is needed, with runs on larger problem sizes but on the same number of resources. In practical terms, we run the B-4 and C-4 instances in addition to the A-4 instance, and thus determine three different instruction rates, one per studied class. Then, depending on whether the current instance handles data that fit in the L2 cache or not, we use the rate coming from the A-4 calibration or the one that corresponds to the instance class.

Note that the use of this cache-aware calibration is mandatory only for the bordereau cluster whose L2 cache is small (1MB per core only). On the graphene cluster that has a L2 cache of 2MB per core, all the instances do fit in cache. Calibrating the simulator with a run of the A-4 instance is then enough.

Impact of the Modifications

Here we evaluate the impact of the modifications detailed in the previous section on each of the issues identified in Section 2.

Instrumentation Time Overhead

Two of the proposed modifications have an impact on the overhead implied by the instrumentation of the application. Applying the -O3 optimization flag transforms the source code and reduces the execution time of both non-instrumented and instrumented versions of the application. Reducing the instrumentation to its minimum will reduce its intrusiveness and thus hinder its impact.

Table 1: Evolution on the bordereau cluster of the execution time and overhead of original and instrumented versions of instances of the LU factorization between the former implementation in [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF] and the modified implementation (-O3 and minimal instrumentation).

Execution times as in [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF] Execution times after modifications Orig.

Instr 1 shows the evolution of the execution times of the original application and its instrumented version between the first implementation presented in [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF] and the current one that implements both modifications. Note that the performance achieved by the current implementation may also be impacted by an update of TAU from version 2.18.3 to version 2.21.1. First, we can note a clear reduction of the execution times of both non instrumented and instrumented versions, coming mainly from the -O3 optimization flag. While simple, this optimization is important as it shortens the overall acquisition time of traces. This gain will become even more interesting for larger problem sizes. If we focus on the evolution of the instrumentation overhead, we see that the impact of the proposed modifications is not stable. This overhead actually increases for two instances (B-8 and C-16) and no clear evolution trend can be found. However, the gain is the most significant for instances that were the most impacted by the instrumentation, justifying the introduction of this modifications. Over the complete set of instances, the instrumentation overhead now ranges from 1.9% to 30% instead of the [4.1%-39.5%] range mentioned in Section 2.

As before, we complete the study on the impact of the proposed modification with results obtained on a more recent compute cluster. We ran experiments with the same initial settings on the graphene cluster, whose results are shown by Table 2.

Table 2: Evolution on the graphene cluster of the execution time and overhead of original and instrumented versions of instances of the LU factorization between the former implementation in [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF] and the current modified implementation (-O3 and minimal instrumentation).

Execution times as in [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF] Execution times after modifications Orig.

Instr As on the bordereau cluster, the use of the -O3 flag has a direct positive impact on the execution time of the application. But the reduction of the instrumentation overhead is more important there. With the first implementation of our framework, the overhead goes from 9.8% to 37.9%, which is consistent with the experiments on the bordereau cluster. The proposed modifications make the overhead range shrink to [0.7% -27.1%]. More interestingly, we can see a stable growth of this overhead as the number of processes increases. As on the bordereau cluster, the improvement is smaller for large number of processes, but remains significant.

Impact of Instrumentation on Hardware Counters Values

To measure the impact of the proposed modifications, i.e., applying the -O3 optimization flag and minimizing the intrusiveness of the instrumentation, on the discrepancy of the measured number of instructions, we made the same experiments as in Section 2.2. Figure 4 shows results obtained on the bordereau cluster. For most instances, the discrepancy is greatly reduced. Except for the B-64 instance, all are now under 6% while they were around 10-13% before. Moreover, we now see the same evolution as in Figure 2 even on the bordereau cluster. This means that the discrepancy now becomes really significant only when the amount of data handled by each process is small. Nevertheless, for the B-64 instance, there is still a significant reduction from 16% to 12% for the worst difference of measurement between both versions. These results are confirmed on the more recent graphene cluster, as shown by Figure 5. The discrepancy is reduced by at least 5% for all instances, most of them being even close to zero. The modifications accentuate the discrepancy increase caused by the evolution of the number of processes. But as for the overhead in terms of execution, instances such as B-128 are not representative of typical executions of Inria parallel applications. Indeed, the benefit of distributing the computation comes to its limits for such instances, with too few data to compute on for each process. Consequently results for these instances are less significant than for instances with a better balance of the computation over communication ratio.

These experiments clearly show the benefits of the proposed modifications on the quality of the instrumentation. We expect it to contribute to the improvement of the accuracy of the simulated execution, studied in the next Section, once combined to a better calibration and the change of simulation back-end.

Accuracy of the Simulated Execution

In this last set of experiments we combine all the proposed modifications. We rely on the compiler optimization flag, the minimal instrumentation, and the cache-aware calibration method to improve the quality of the simulation of the compute part of the application. For the communication part, we benefit of the quality of the underlying SMPI layer that provides good estimations for the point-to-point communication of small messages and a tuned piece-wise linear network model.

Figure 6 shows the impacts of the whole set of modifications on the accuracy of the simulated executions of the LU factorization. These results, obtained on the bordereau cluster, have to be compared to those presented by Figure 3. The first outcome is a drastic reduction of the inaccuracy. With our initial implementation, the relative error linearly increased from -2.7% to 38.9% for Class B and from -15.8% to 32.5% for Class C. Now, the error is comprised between -9.5% and 11.5% for Class B and between -6.1% and 8.1% for Class C. Moreover, the linear increase of the error along with the number of processes, that indicated a bad simulation of the communication part, disappeared. We even note an opposite trend for Class B. This can be explained by the fact that SMPI does not model the time to copy data in memory in the MPI_Send function yet. This leads to a slight underestimation of the simulation time that increases as more small messages are exchanged for large number of processes. This phenomenon is well identified and will soon be addressed.

Figure 7 shows the relative error achieved under the same experimental conditions on the graphene cluster. As on the bordereau cluster, the error is in a narrow interval ranging from -11.4% and -2%. Results for Class B are even more stable. Again the underestimation of the simulated time should be compensated by taking memory copy into account.

While narrowing the inaccuracy range is a good result in itself, the more important result shown by Figures 6 and7 is the greater stability of this relative error. Indeed, the strong linear increase highlighted by Figure 3 prevented us to use our framework for its main objective, i.e., predicting the execution time of any instance of a given application within a certain confidence interval.

The results presented in this section are not perfect yet, as some fluctuations or linear trends are still there, but there are a great improvement towards an accurate and trustworthy performance prediction tool. Some issues are still to be addressed, but most of them have been identified. Finally the factorization of the code, and thus of the efforts with those made on the SMPI project will help to improve the quality of our Time-Independent Trace replay tool.

Related Work

The off-line simulation of MPI application has been used extensively, as shown by the number of tracebased simulators described in the literature since 2009 [START_REF] Hoefler | LogGOPSim -Simulating Large-Scale Applications in the LogGOPS Model[END_REF][START_REF] Tikir | PSINS: An Open Source Event Tracer and Execution Simulator for MPI Applications[END_REF][START_REF] Núñez | New Techniques for Simulating High Performance MPI Applications on Large Storage Networks[END_REF][START_REF] Zhai | PHANTOM: Predicting Performance of Parallel Applications on Large-Scale Parallel Machines Using a Single Node[END_REF][START_REF] Hermanns | Verifying Causality between Distant Performance Phenomena in Large-Scale MPI Applications[END_REF].

The typical approach is to decompose the trace in time intervals delimited by the MPI communication operations. The application is thus seen as a succession of computation and communications steps. Most of the existing tools then simply replace the computations by the measured delays. The performance differential between the platform used to obtain the trace and the target platform can be taken into account, but usually using a simple scaling [START_REF] Tikir | PSINS: An Open Source Event Tracer and Execution Simulator for MPI Applications[END_REF][START_REF] Núñez | New Techniques for Simulating High Performance MPI Applications on Large Storage Networks[END_REF][START_REF] Hermanns | Verifying Causality between Distant Performance Phenomena in Large-Scale MPI Applications[END_REF]. Network communications are simulated based on the communication events recorded in the trace and on a simulation model of the network. As for computation, most of the tools adopt simplistic network models. The most common simplifications are to ignore network contention because it is known to be costly to simulate [START_REF] Zheng | Simulation-Based Performance Prediction for Large Parallel Machines[END_REF], and use monolithic performance models of collective communications rather than simulating them as sets of point-to-point communications [START_REF] Tikir | PSINS: An Open Source Event Tracer and Execution Simulator for MPI Applications[END_REF][START_REF] Badia | Dimemas: Predicting MPI applications behavior in Grid environments[END_REF]. Exception is the MPI-NetSim on-line simulator [START_REF] Penoff | MPI-NeTSim: A network simulation module for MPI[END_REF], which relies on a slow packet-level discrete event network simulator.

Inria 6 Conclusion and Future Work

Simulation is a popular approach to obtain objective performance indicators platforms that are not at one's disposal. It may help the dimensioning of compute clusters in large computing centers. In [START_REF] Desprez | Assessing the Performance of MPI Applications Through Time-Independent Trace Replay[END_REF], we proposed a new approach for the off-line simulation of MPI applications. Instead of relying on logs of execution that associate an event to a time-stamp, we use Time-Independent Traces as an input of our simulator. These traces contain only information about volumes of computation and communications. This allows us to decouple the acquisition process from the replay of the trace. Heterogeneous and distributed platforms can then be used to get traces without impacting the quality of the simulation, which is not possible with any other tool. Our simulator is built on top of the SIMGrid toolkit and thus benefits of its fast, scalable, and validated simulation kernel. We also rely on a well established tracing tool, TAU, in the acquisition process.

In this paper, we identified several issues, i.e., an instrumentation time overhead, a discrepancy of the measured number of instructions, a calibration method that ignores some cache related effects, and an improvable simulation back-end, related to the first implementation of this tool. These issues lead to inaccurate and unstable simulated times and prevent a use for performance prediction within a confidence interval. We proposed several modifications of our framework to address these different issues and evaluated their respective impact. While not perfect, the resulting implementation is far more accurate and stable than the previous one, and becomes a good candidate to predict the performance of MPI applications on platforms that are not available.

This Time-Independent Trace replay framework is freely distributed as part the SIMGrid project under the LGPL license. It can be downloaded from http://simgrid.gforge.inria.fr/ download.php. The implementation described in this article will be released in the version 3.8 of SIMGrid. Moreover, the process to acquire a Time-Independent Trace has been fully described in [START_REF] Markomanolis | Time-Independent Trace Acquisition Framework -A Grid'5000 How-to[END_REF].

As future work, we plan to implement the missing feature to model the time taken in sends and receives to copy data in memory in the eager mode of MPI. We also aim at improving our calibration method to automatically take cache usage into account and better estimate the instruction rate used by the simulator. Finally, we would like to assess the performance of the proposed approach on an MPI application used in production.

Figure 1 :

 1 Figure 1: Distribution across processes of the relative difference (in %) of the measured numbers of executed instructions between fine and coarse grain instrumented versions of various instances of the LU factorization on the bordereau cluster.

Figure 2 :

 2 Figure 2: Distribution across processes of the relative difference (in %) of the measured numbers of executed instructions between fine and coarse grain instrumented versions of various instances of the LU factorization on the graphene cluster.

  time (in %) Number of processes

Figure 3 :

 3 Figure 3: Evolution with regard to the number of processes of the relative error between execution and simulated times for the execution of the LU factorization.

Figure 4 :

 4 Figure 4: Distribution across processes of the relative difference of the measured numbers of executed instructions between minimal and coarse instrumentations of optimized (-O3) LU instances. Results obtained on the bordereau cluster.

Figure 5 :

 5 Figure 5: Distribution across processes of the relative difference of the measured numbers of executed instructions between minimal and coarse instrumentations of optimized (-O3) LU instances. Results obtained on the graphene cluster.

Figure 6 :

 6 Figure 6: Evolution w.r.t. the number of processes of the relative error between execution and simulated times for the execution of the LU factorization with the new replay framework. Results obtained on the bordereau cluster.

Figure 7 :

 7 Figure 7: Evolution w.r.t. the number of processes of the relative error between execution and simulated times for the execution of the LU factorization with the new replay framework. Results obtained on the graphene cluster.

RR n°8092

Acknowledgments

This work is partially supported by the ANR project SONGS (11 ANR INFRA 13) and the CNRS PICS N • 5473. Experiments presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RENATER and several Universities as well as other funding bodies (see https://www.grid5000.fr). The authors would also like to thanks all the members of the SIMGrid project working on the SMPI effort, and especially A. Degomme and L. Schnorr for their great help.


Benjamin Marty

MSc Benoit Larrat

Maxime Van Landeghem

PhD Caroline Robic

PhD Philippe Robert

PhD Marc Port

PhD Denis Le Bihan

MD, PhD Mathieu Pernot

PhD Mickael Tanter

PhD Franck Lethimonnier

PhD Maxime Van Landeghem

PhD Sébastien Mériaux
email: sebastien.meriaux@cea.fr

Dynamic study of blood-brain barrier closure after its disruption using

Keywords: MRI, Molecular imaging, Pharmacology, Ultrasound, Blood brain barrier BBB disruption, MR contrast agent, MRI, T 1 mapping, ultrasound, nanoparticles

Introduction

Cerebral tissues are isolated from circulating blood by the blood brain barrier (BBB) [START_REF] Rubin | The cell biology of the blood-brain barrier[END_REF]. This physiological barrier consists of a lining of tightly packed vascular endothelial cells, different from the peripheral blood vessels [START_REF] Hawkins | The blood-brain barrier/neurovascular unit in health and disease[END_REF]. The tight junctions between these cells restrict the diffusion of microscopic objects (e.g. bacteria) and of large, hydrophilic molecules (>400 Da) from blood to brain parenchyma, while allowing the passage of small, hydrophobic molecules such as O 2 , CO 2 , proteins, and metabolites [START_REF] Pardridge | The blood-brain barrier: bottleneck in brain drug development[END_REF]. Because of the neuroprotective nature of the BBB, delivery of potentially important diagnostic and therapeutic agents is a major challenge in the treatment of most brain disorders. Strategies to design specific drugs targeted to the brain involve finding an efficient mode of drug delivery across the BBB.

Recently, it has been demonstrated that the use of low power focused ultrasound combined with a systemic injection of lipid-(or polymer-) shelled microbubbles enables a non-invasive, local and transient disruption of the BBB [START_REF] Hynynen | Noninvasive MR imagingguided focal opening of the blood-brain barrier in rabbits[END_REF]. Many studies were then carried out to (i) establish optimal ultrasound parameters that permit adequate tissular penetration without causing tissue damage [START_REF] Choi | Optimization of blood-brain barrier opening in mice using focused ultrasound[END_REF][START_REF] Hynynen | Local and reversible blood-brain barrier disruption by noninvasive focused ultrasound at frequencies suitable for trans-skull sonications[END_REF][START_REF] O'reilly | The impact of standing wave effects on transcranial focused ultrasound disruption of the blood-brain barrier in a rat model[END_REF]O'Reilly et al 2011a;O'Reilly et al 2011b;Sheikov et al 2008), (ii) quantify permeability of the disrupted brain tissue [START_REF] Vlachos | Permeability assessment of the focused ultrasound-induced blood-brain barrier opening using dynamic contrast-enhanced MRI[END_REF][START_REF] Vlachos | Permeability Dependence Study of the Focused Ultrasound-Induced Blood-Brain Barrier Opening at Distinct Pressures and Microbubble Diameters Using DCE-MRI[END_REF], and (iii) evaluate responses to treatments of particular brain disorders including tumors [START_REF] Chen | Novel magnetic/ultrasound focusing system enhances nanoparticle drug delivery for glioma treatment[END_REF][START_REF] Liu | Blood-Brain Barrier Disruption with Focused Ultrasound Enhances Delivery of Chemotherapeutic Drugs for Glioblastoma Treatment[END_REF][START_REF] Treat | Impact of Focused Ultrasoundenhanced Drug Delivery on Survival in Rats with Glioma[END_REF] and Alzheimer's disease [START_REF] Jordao | Antibodies targeted to the brain with image-guided focused ultrasound reduces amyloid-beta plaque load in the TgCRND8 mouse model of Alzheimer's disease[END_REF][START_REF] Raymond | Ultrasound Enhanced Delivery of Molecular Imaging and Therapeutic Agents in Alzheimer's Disease Mouse Models[END_REF]. Most of these studies used magnetic resonance contrast agents (MR-CA) for monitoring the processes. Despite a rapidly growing number of studies, the mechanism of ultrasound-induced BBB opening is understood only poorly. In particular, the maximum space that can be safely generated (ensuring reversibility) between endothelial cells, and the duration for which this 4 opening lasts, have not been measured. The knowledge of these parameters is crucial for the current development of brain targeted nanoparticles. Functionalized MR-CA for instance, are available in a wide range of hydrodynamic diameters: gadolinium chelates, a few nm [START_REF] Wadghiri | Detection of Alzheimer's amyloid in Transgenic mice using magnetic resonance microimaging[END_REF], iron nanoparticles, 25-100 nm (Mendonca [START_REF] Dias | Ferromagnetic particles as contrast agents for magnetic resonance imaging of liver and spleen[END_REF][START_REF] Renshaw | Ferromagnetic contrast agents -A new approach[END_REF], gadolinium based emulsions or liposomes, 200-300 nm [START_REF] Devoisselle | Entrapment of Gadolinium-DTPA in liposomes -Characterization of vesicles by P-31 NMR-Spectroscopy[END_REF], and have very different vascular remanences, ranging from a few minutes to tens of hours. An estimate of the time window during which a molecule of a given size may be delivered across the BBB would be useful to monitor the amount of drug released to cerebral tissue and to adjust the dosage. For example, it has been shown that 3 -8 nm wide particles are able to cross BBB 20 min after its ultrasound-induced disruption, but not particles of size 50 nm [START_REF] Choi | Molecules of various pharmacologically-relevant sizes can cross the ultrasound-induced Blood-Brain Barrier opening in vivo[END_REF]. Additionally, as the brain is unprotected while the BBB is open, it is important to know the time-to-recovery to limit any possible brain damages induced by tissular penetration of pathogenic agents.

The quantification of the amount of MR-CA crossing the BBB requires development of molecular imaging sequences via imaging techniques that are sensitive (high field MRI) and provide outputs that can be correlated with the MR-CA concentration. Contrast agent (CA) quantification is one of the principal challenges for MR molecular imaging. Unlike other techniques such as fluorescence or nuclear imaging where the received signal comes directly from the CA, MR signal comes from the surrounding water protons and is therefore indirectly related to the presence of an exogenous probe. To quantify the CA concentration using MRI, it is necessary to model its interaction with the surrounding water molecules.

Paramagnetic complexes interact with the surrounding water protons by decreasing their relaxation times [START_REF] Swift | NMR-Relaxation mechanisms of 017 in aqueous solutions of paramagnetic cations and lifetime of water molecules in first coordination sphere[END_REF]. CA concentration can be linked to the T 1 or T 2 relaxation time decrease (in ms). In this study, we adapted a T 1 mapping sequence proposed by Deichmann and colleagues [START_REF] Deichmann | Quantification of T1 values by Snapshot-Flash NMR imaging[END_REF][START_REF] Deichmann | Fast T1 mapping on a whole-body scanner[END_REF] In this study, a calibrated BBB opening procedure was performed in healthy young adult rats under MR guidance. Five MR-CA were used (3 paramagnetic and 2 superparamagnetic) with different hydrodynamic diameters (from 1 to 65 nm) to investigate the maximum permissible size of the gaps induced in the endothelial wall under safe conditions.

Then, using the T 1 mapping strategy for the 3 paramagnetic CA we quantified the amount of particles crossing the BBB when injected at different times after ultrasound-induced BBB disruption. This enabled quantitative monitoring of the dynamics of BBB closure for each given molecular size. Lastly, we present a theoretical model to fit the experimental data and derive a calibration curve to predict duration of BBB opening as a function of the hydrodynamic diameter of a given CA.

Materials and Methods

Animal Preparation

All experiments were performed in accordance with the recommendations of the European Community (86/609/EEC) and the French legislation (decree no. 87/848) for use and care of laboratory animals. A total of 36 Sprague Dawley male rats (80-100 g, Janvier, France) were used. Their head was shaved to ensure proper coupling of the ultrasound beam to the brain. Rats were anesthetized with 1.5% isoflurane in a mixture of air and oxygen and placed in a cradle in prone position. A catheter (25 G needle) was positioned in the caudal vein to inject microbubbles and MR-CA from outside the scanner with minimal dead volume. Table 1. After each MR-CA injection, the bolus was flushed by injection of 100 µl of saline solution. Three of the five MR-CA were paramagnetic Gd-chelates (Dotarem®, P846, and P792) and were detected and quantified using T 1 strategies [START_REF] Kang | Non-Invasive Visualization of Basilar Artery Perforators With 7T MR Angiography[END_REF][START_REF] Protti | Late Gadolinium Enhancement of Acute Myocardial Infarction in Mice at 7T: Cine-FLASH Versus Inversion Recovery[END_REF][START_REF] Yankeelov | Repeatability of a reference region model for analysis of murine DCE-MRI data at 7T[END_REF]. The other two CA (P904 and P03680) were USPIO. A T 2 -weighted sequence was use to detect them in cerebral tissues [START_REF] Hyodo | Estimation of tumor microvessel density by MRI using a blood pool contrast agent[END_REF][START_REF] Philippens | Characterization of late radiation effects in the rat thoracolumbar spinal cord by MR imaging using USPIO[END_REF].

Longitudinal r 1 and transverse r 2 relaxivities were measured in vitro at 7T using galleries of tubes containing different concentrations of CA diluted in a 0.3% agar matrix, and maintained at 37°C.

Ultrasound equipment

A MR-compatible focalized transducer (central frequency 1.5 MHz, diameter 30 mm, focal depth 20 mm, Imasonic, France) driven by a programmable function generator was used to produce ultrasound waves. It was coupled to the rat skull using a latex balloon filled with deionized and degassed water. Electrical power sent to the transducer was monitored during the BBB opening session. Main characteristics of the transducer (focal point size 0.6×0.6×3 mm 3 , transcranial acoustic transmission factor 52 ± 5 %) were estimated in a previous study [40].

BBB opening

BBB disruption was performed during MR imaging session thanks to a dedicated holder maintaining the ultrasound transducer above rat head. Its position was monitored using the ARFI sequence and the right thalamus in the brain was chosen as a target for BBB bubbles/ml, 200 µl, 2 s) approximately 5 s before the beginning of ultrasound session. The bolus was flushed by injection of 100 µl of saline solution. Sonication was performed with 3 ms bursts every 100 ms for one minute [START_REF] Choi | Noninvasive and localized blood-brain barrier disruption using focused ultrasound can be achieved at short pulse lengths and low pulse repetition frequencies[END_REF]. The peak negative acoustic pressure at the ultrasound focal point was calibrated to 0.45 MPa, a level sufficiently low to ensure safe and reversible BBB opening [START_REF] Larrat | MR-Acoustic Radiation Force Mapping can Quantitatively Predict Drug Delivery following Ultrasound-Induced Blood Brain Barrier Disruption in Rats at High Field)[END_REF][START_REF] Chopra | Influence of exposure time and pressure amplitude on bloodbrain-barrier opening using transcranial ultrasound exposures)[END_REF]. Anatomical images were acquired at the end of each MRI session to verify integrity of the brain tissues.

MRI acquisitions

MRI was performed on a 7T/90 mm Pharmascan scanner (Bruker, Germany). A saddle coil was specially designed in-house for excitation and signal reception. The geometry and size of the resonator (diameter 30 mm) were optimized to maximize the overall signal-tonoise ratio over the whole brain while allowing the ultrasound beam to propagate from the transducer to the rat brain.

A multislice spin-echo sequence was modified to include additional motion-sensitizing gradients, and synchronized to ultrasonic bursts so that the phase signal was proportional to local acoustic intensity [START_REF] Larrat | MR-guided transcranial brain HIFU in small animal models[END_REF][START_REF] Mcdannold | Magnetic resonance acoustic radiation force imaging[END_REF]. This ARFI sequence was acquired with the following parameters: TE/TR = 40/1700 ms, T acq = 4 min, R = 0.5×0.5×1 mm 3 , duration of the motion encoding gradients = 13.3 ms, duration of the sonication = 3 ms.

A high spatial resolution T 2 -weighted RARE sequence (TE eff /TR = 32/4200 ms, R = 0.125×0.125×0.3 mm 3 ) was used to detect the presence of USPIOs, and also for acquiring images at the end of every experiment to confirm lack of haemorrhages or edema due to To measure concentration of Gd-chelates, a T 1 mapping sequence was acquired before and after MR-CA injection. It consisted in a segmented series of fast gradient echo (FGE) images acquired at different time points after magnetization inversion in order to follow the entire T 1 recovery curve [24,25]. Sequence parameters were: TR 1 , 5 ms, TE, 2.5 s, 6 segments, 60 inversion times (from 64 to 5800 ms), flip angle, 5°, and R = 0.2×0.2×1 mm 3 .

Repetition time between the acquisitions of two segments, TR 2 , was 9 s, and total acquisition time is 12.5 min. A centric encoding of the k-space was chosen in order to avoid T 2 effects during echo trains at the acquisition of central lines.

Data analysis

Data were analyzed using dedicated codes written in Matlab software (MathWorks, USA). Maps of acoustic pressure index (IP ac ) were deduced from ARFI sequence using the following equation:

th OFF ON ac P IP ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ - -- - ∝ ∝ ∝ ∝ (1)
where φ ON (respectively φ OFF ) is the phase of ARFI signal acquired with ultrasound (respectively without ultrasound), and P th is the expected acoustic pressure at the focal point in case of an ideal ultrasound coupling. In a previous study [START_REF] Larrat | MR-guided transcranial brain HIFU in small animal models[END_REF], we verified the linearity between acoustic intensity and MR-ARFI phase signal using the same setup. In another study [START_REF] Larrat | MR-Acoustic Radiation Force Mapping can Quantitatively Predict Drug Delivery following Ultrasound-Induced Blood Brain Barrier Disruption in Rats at High Field)[END_REF], we also characterized and verified the linearity between acoustic pressure and the amount of MR-CA crossing the BBB in the range of acoustic pressure values applied here. To generate T 1 maps, the MRI signal measured with the FGE sequence was fitted as a function of the inversion time pixel by pixel as proposed by Deichmann and colleagues [24,25]. MR-CA concentration maps (C) were then calculated from the T 10 and the T 1 maps after MR-CA injection using the following equation [START_REF] Swift | NMR-Relaxation mechanisms of 017 in aqueous solutions of paramagnetic cations and lifetime of water molecules in first coordination sphere[END_REF], considering that relaxivities r 1 measured in vitro in agar matrix hardly differ from the ones in rat brain tissues:

                                - -- - = = = = 10 1 1 1 1 1 T T r C
( In a rat model, BBB was opened transiently with the use of ultrasound, immediately followed by intravenous injection of Dotarem® (gadolinium (Gd)-chelate) and measurement of the amount of the MR-CA delivered across the BBB (Figure 1). Before the ultrasound sonication, acoustic radiation force imaging (ARFI) was performed to obtain the acoustic pressure map in the brain (Figure 1-A). This sequence was used to verify that the ultrasound beam was focused at the desired point in the right thalamus, and to estimate the acoustic pressure at this focal point prior to microbubbles injection, ensuring that the acoustic pressure remains below a predetermined safety threshold [START_REF] Chopra | Influence of exposure time and pressure amplitude on bloodbrain-barrier opening using transcranial ultrasound exposures)[END_REF][START_REF] Larrat | MR-Acoustic Radiation Force Mapping can Quantitatively Predict Drug Delivery following Ultrasound-Induced Blood Brain Barrier Disruption in Rats at High Field)[END_REF]. Once the transducer was set at the desired position, Sonovue® microbubbles were injected intravenously via a catheter positioned in the caudal vein, and the BBB was disrupted by pulsed sonication for one minute. Dotarem® was then injected using the same catheter. To measure the concentration of Dotarem® delivered, T 1 maps were acquired before and after MR-CA injection (Figure 1-B). There was high tissue perfusion of the CA (strong T 1 decrease) in the cheek and neck muscles, but not in the brain tissue due to the BBB. As expected, the ultrasound focal point (red arrow) was the only area of the brain exhibiting a strong T 1 decrease indicating BBB disruption. Based on T 1 mapping, CA concentration maps were obtained at the three selected regions of interest (ROIs): ultrasound focal point, ROI spot (n°1), a control region contralateral to the ultrasound focal point, ROI contra (n°2) and a region in cheek muscles ROI cheeks (n°3) (Figure 1-C).

Maximum gap between endothelial cells after reversible ultrasound-induced disruption of the BBB

In order to determine the size of gaps generated after ultrasound-induced disruption of the BBB, we chose five different MR-CA of varying hydrodynamic diameters ( 11 injected (n=1 rat for each CA) directly after ultrasound-induced BBB disruption. T 1 -weighted (T 1w ) images were obtained for the 3 paramagnetic CA and T 2 -weighted (T 2w ) images for the two ultrasmall superparamagnetic iron oxide (USPIO) CA (Figure 2). For all three Gdchelates, a positive contrast was clearly visible on T 1w images at the ultrasound focal point, revealing the presence of the CA. Student's t-test performed on the measured MR signals confirmed that signal differences between ROI spot and ROI contra were significant (p < 0.01).

For the 25 nm USPIO, a negative contrast appeared on T 2w images after its injection suggesting that BBB was permeable to molecules of this size. Signal differences with contralateral region were also significant (p < 0.01). However, the largest USPIO (65 nm) was detected only at the very focal point of the transducer as its penetration was significantly hindered compared to the smaller USPIO. The signal differences between ROI spot and ROI contra for this USPIO were not significant (p > 0.01) if the ROI was of the standard size (0.6×0.6×3 mm 3 ), but the differences became significant (p < 0.01) if the size of the ROI spot was decreased to 0.2×0.2×1 mm 3 around the focal point. These results suggest that 65 nm is close to the maximal gap between endothelial cells that is attainable after BBB disruption with this protocol.

Quantification of MR-CA concentration

Rats were injected with Dotarem® at increasing doses (375, 750 and 1500 µmol/kg, n=1 rat per dose) following BBB disruption. Concentration maps acquired by T 1 mapping were overlaid on the anatomical T 2w images (Figure 3-A). Additionally, a T 1 -weighted MSME sequence was also acquired. The estimated CA concentration at the ultrasound focal point was found to be proportional to the injected dose (C 375 < C 750 < C 1500 ) (Figure 3-B). In contrast, the T 1w signals did not exhibit the same behavior and even saturated at high CA injected doses (Figure 3-C). These observations confirm the high interest of our protocol to

C o n f i d e n t i a l : F o r R e v i e w O n l y 12
precisely quantify small variations of CA concentration, and that quantitative analysis based on T 1w signals alone can be misleading. Hence, results presented in the following sections of this study were all derived from quantitative data obtained using the T 1 mapping sequence (Equation 3).

Quantitative study of BBB closure dynamics

The BBB closure dynamics was assessed using paramagnetic CA of different molecular sizes, administrated at different times after ultrasound-induced disruption. Rats were injected with the following CA: Dotarem® (n=10 rats), P846 (n=10 rats) and P792 (n=8 rats) from 0 to 24 hours after ultrasound sonication. The reproducibility of CA delivery was verified by measuring both the injected dose and the acoustic power deposited at the ultrasound focal point. CA concentrations measured in the cheek muscles in each animal were directly correlated with the injected dose since no permeability barrier like BBB exists in muscles (Figure 4-A). Mean values were: 0.345 ± 0.034 mM for Dotarem® group, 0.061 ± 0.008 mM for P846 group and 0.040 ± 0.002 mM for P792 group. Variabilities observed within the above three groups (10%, 13% and 5%, respectively) could be explained by experimental errors during injection and also by slight variations in the weight or metabolism of the animals. An index of the acoustic pressure (IP ac ) deposited at the focal point was estimated from ARFI sequence (Figure 4-B) and was 1.13 ± 0.06 for Dotarem® group, 1.00 ± 0.05 for P846 group and 0.97 ± 0.13 for P792 group. Variabilities observed within the above three groups (5%, 5% and 13%, respectively) was mainly explained by slight differences in the placement of the transducer on the heads. Variations in the incidence angle of the ultrasound beam or in the thickness of the skull could have a significant impact on attenuation of the beam during its passage through the skull and lead to variations of the peaknegative acoustic pressure at the focal point. For Dotarem® (1 nm particle), the BBB was permeable for several hours post-BBB disruption. A significant amount of CA was detectable in the brain even when injected 8 or 9 hours after disruption. Around 24 h post-disruption, the BBB had apparently recovered its whole integrity, as very few CA were detectable in brain tissues. P846 (4 nm particle) was delivered across the BBB only if injected within 2 hours post-BBB disruption. For P792 (7 nm particle), the duration of CA delivery to the brain was even shorter than for P846. One hour post-BBB disruption, the amount of CA crossing the BBB was about 30% of its initial value and dropped to undetectable levels after 90 minutes.

These findings confirmed that the duration of BBB crossing decreases when increasing CA hydrodynamic diameter. Secondly, it appears from the shape of the timedependent concentration curve (Figure 4-C) that the BBB permeability decreased at a faster rate initially and at a lesser rate later on.

A model for BBB closure dynamics

Based on the experimental observations, the BBB closure dynamics was modelled as follows. First, it was assumed that each individual gap generated in the BBB by ultrasound returns to its equilibrium position (i.

- -- - - -- - - -- - ∝ ∝ ∝ ∝ σ σ σ σ σ σ σ σ (4)
For each time t, the amount Q CA (t) of CA of hydrodynamic diameter d H crossing the BBB in a unit volume post-disruption is proportional to the inflow flux only, since the outflow of CA from brain to blood is assumed to be negligible due to high concentration in the blood compartment. Assuming a hard sphere behavior for CA molecules, Q CA (t) can be expressed as follows:

∫ ∫ ∫ ∫ +∞ +∞ +∞ +∞ ∝ ∝ ∝ ∝ H d H CA dx d x t x N t Q 2 ). , ( ) ( (5)
The CA diffusion in brain outside of the BBB disrupted region is also neglected since the CA tissular concentration is quantified within a few minutes after CA intravenous injection. A previous study from our team [START_REF] Marty | Apparent diffusion coefficient of Gd-based contrast agents assessed in vivo in the rat brain using dynamic T1 mapping)[END_REF] showed that the Gd-chelates used in this study do not diffuse fast enough in brain parenchyma to make this effect significant: for example, the diffusion coefficient of Dotarem® (1 nm particle) was estimated around 46 µm 2 .s -1 (Table 1). The half closure time (t 1/2 ) was then defined as the time after disruption when the concentration of CA crossing the BBB was 50% of the maximal concentration obtained immediately after disruption. t 1/2 was then calculated as the numerical solution of the following equation:

Integration of equation 5 gives:

                                        + + + +                                                                 - -- - ∝ ∝ ∝ ∝ - -- - - -- - - -- - - -- - - -- -
C o n f i d e n t i a l : F o r R e v i e w O n l y 15                                         + + + +                                                                 - -- - = = = =                                         + + + +                                                                 - -- - - -- - - -- - - -- - - -- - - -- - - -- - 2 / 1 2 2 0 2 2 / 1 2 / 1 2 / 1 2 0 2 2 0 0 2 2 0 0 2 1 2 2 1
σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ σ π π π π σ σ σ σ σ σ σ σ π π π π (7)
From the fitted parameters (σ 0 and k), half closure times could be estimated as a function of hydrodynamic diameters of the molecules by numerically solving equation 7 (Figure 5, black circles). The calculated half closure time was ~ 5.5 h for a 1 nm particle; ~ 1.5 h for a 4 nm particle and about 30 min for a 7 nm particle. These values were close to those estimated experimentally for Dotarem®, P846 and P792, which were 5.2 h, 1.2 h and 45 min, respectively (Figure 5, red crosses).

For larger molecules, half closure times were greatly reduced (around 3 min for a 25 nm particle and only few seconds for a 65 nm particle). These results match the qualitative results obtained with USPIOs (Figure 2). Therefore, to deliver a significant amount of larger molecular size CA across the BBB their intravenous circulation should coincide with the end of ultrasound sonication.

Lastly, it appears that the following analytical formula fits reasonably well the numerical solution of equation 7, enabling to directly estimate half closure time as a function of hydrodynamic diameter:

2 2 / 1 . 1 H d B A t + + + + = = = = (8)
with A = 2.35×10 4 s and B = 0.2106 nm -2 (Figure 5, solid line).

Discussion

Ultrasound-induced BBB disruption

Drug delivery to the brain remains a major challenge for the treatment of most brain diseases. Although pathologies such as cancer degrade BBB integrity (Liu et

MR guidance

The ultrasound-induced BBB disruption described here was performed entirely under high-field MRI guidance. This is a challenge as it requires integrating a stereotactic frame, a MR-compatible ultrasound transducer adapted to the geometry of the rodent head and a dedicated radiofrequency (RF) coil fitted into the 8.5 cm bore of a preclinical MRI scanner.

The optimized RF coil allowed the quantification of Gd MR-CA at a micromolar range, with sufficient signal homogeneity over the whole brain and the ability for the ultrasound beam to propagate through the coil.

The anesthetized animal was placed inside the MRI scanner just before the start of the BBB disruption procedure till recovery. Thus, it was possible to obtain an acoustic intensity map before BBB disruption. 17 irreversible tissue damage. Furthermore, microbubbles were injected to cause ultrasoundinduced BBB disruption followed by injection of MR-CA via the same the catheter, which allowed imaging of the early events post-BBB disruption.

Maximum gap obtained between endothelial cells

CA of increasing hydrodynamic sizes were tested, and the maximum gap width generated between endothelial cells after BBB disruption was estimated with our acoustic parameters. Gaps slightly above 65 nm could be generated safely, and the BBB recovered full integrity within 24 h. The mechanical index of our set-up was 0.37 which is below the standard value proposed by other studies for reversible BBB opening (McDannold et al 2008). This result is particularly interesting from the view of development of targeted CA.

Using a similar sonoporation protocol, it would be possible to deliver nanoparticles such as USPIOs, functionalized to target tissular biomarkers. In contrast, larger objects like emulsions or liposomes that typically have hydrodynamic diameters >100 nm, would be difficult to deliver to the brain using this technique, assuming a purely paracellular pathway. In general, an estimate of the maximum possible gap diameter would be helpful when the technique is applied to deliver therapeutic molecules. Several therapeutic drugs could benefit from this information: for instance chemotherapeutic drugs such as Avastin® An optimization of parameters based on the knowledge gained from these studies may allow generation of pores larger than 65 nm wide.

BBB closure dynamics

In this study, the concentration of paramagnetic CA delivered at the ultrasound focal point was quantitatively measured, based on which, the duration of passage across the BBB was estimated for molecules of different sizes. This duration decreased rapidly as a function of hydrodynamic diameter of the CA. For instance, small molecules (around 1 nm) continued to cross the BBB for more than 10 h, whereas large iron oxide particles (around 25 nm) had a time window of only a few minutes. Therefore, to maximize the amount of extravasation in the brain, large molecules should be injected right after or even during ultrasound sonication.

The results presented here also suggest a closure of the BBB at a progressively decreasing rate. We present a simple, damped, elastic model assuming an exponential decay of individual gap sizes to semi-quantitatively describe the observed dynamics of BBB closure.

The experimental data were fed in the model to derive important physiological parameters, such as the typical closure time constant k that was estimated to be 1.54×10 -5 s -1 . This parameter is an indirect measurement of rate of endothelial cell relaxation after contraction. It could be used as a physio-pathological index of vessel integrity.

This model assumes a paracellular passage of injected molecules with an inflow limited by the size of the molecule relative to the size of the pores. Our results support the notion that most of the crossing happens between endothelial cells rather than through them.

This reinforces the hypothesis that during ultrasound-induced BBB opening sonicated

C o n f i d e n t i a l : F o r R e v i e w O n l y
19 microbubbles exert a mechanical stress on endothelial cells which makes them contract on themselves and loosen the junctions between them creating gaps [START_REF] Vykhodtseva | Progress and problems in the application of focused ultrasound for blood-brain barrier disruption[END_REF].

Limitations of the study

The proposed BBB closure model has certain limitations. It does not take into account potential increased permeability of endothelial cell membranes (transcellular BBB crossing).

It does not distinguish between hydrophilic and lipophilic molecules although it is known that lipophilic agents present an increased ability to cross the BBB via transcytosis [START_REF] Pardridge | The blood-brain barrier: bottleneck in brain drug development[END_REF]. This concurrent release pathway could explain the observed differences between experimental data of P846 and P792 and predictions from the theoretical model.

The three Gd-chelates compared in this study have similar but not equal vascular remanences (Table 1). The tissular concentrations were measured within 12.5 minutes after CA injection, which is not a negligible time delay as the plasmatic half-life of the CA are between 26 and 51 min [START_REF] Fries | Brain Tumor Enhancement in Magnetic Resonance Imaging at 3 Tesla Intraindividual Comparison of Two High Relaxivity Macromolecular Contrast Media With a Standard Extracellular Gd-Chelate in a Rat Brain Tumor Model[END_REF]. Although during data processing the differences in overall injected plasmatic concentration were corrected for each CA (Equation 3), temporal variations of vascular concentrations during the T 1 mapping acquisition could be different from one CA to the other, none of them remaining strictly steady. The impact of these variations may be low since the MRI mapping acquisition was performed with a centric encoding scheme. Nevertheless, it may explain the differences between the theoretical model and the experimentally measured closure dynamics for the different molecules.

Conclusion

In this study, we developed a complete methodology for noninvasive, reversible and controlled ultrasound-induced BBB disruption under MR-guidance that permitted precise Based on the theoretical model proposed here, half closure time, t 7) (black circles). Red crosses correspond to experimental data on the paramagnetic nanoparticle contrast agents used here (Dotarem® 1 nm; P846 4 nm; and P792 7 nm) and the solid line represents the best fit with analytical function: t1/2 = A / (1 + B.dH 2 ), where dH is the hydrodynamic diameter of the injected CA (in nm) and A and B are constants, A = 2.35×10 4 s and B = 0.2106 nm -2 . 80x68mm (300 x 300 DPI)

  10% heparin was added to all injected solutions to avoid clot formation in the catheter. Body temperature and respiration rate were continuously monitored during the experiments. of different hydrodynamic diameters were provided by Guerbet Research (France). Main characteristics and injected doses of these nanoparticles are summarized in

  microbubbles (Bracco, Italy) were administrated via a bolus (1.5 × 10 8

  weighted MSME sequence (TE/TR = 8/300 ms, R = 0.250×0.250×1 mm 3 ) was acquired to detect paramagnetic Gd-chelates.

  2) Signals from different manually drawn ROIs were analyzed. A ROI spot of 0.6×0.6×3 mm 3 corresponding to the transducer focal point size was drawn in the right thalamus. To analyze the effect of BBB disruption, a similar ROI contra was taken in the corresponding contralateral (left thalamus) region of the brain. This enabled correction for the residual vascular concentration of MR-CA during imaging. Another control ROI was drawn in an extracerebral region in the cheek muscles (ROI muscle ). For each rat, a corrected MR-CA concentration at the focal point (C* without unit) was calculated as given below, based on the concentration measured at the ROI spot (C ROIspot ) and that at the control spots (C ROIcontra and C ROImuscle ), and the index of deposited acoustic pressure (IP ac ): was mandatory to ensure the correction of the bias introduced by unequal injected volumes (unequal plasmatic concentrations) and the bias introduced by variations in acoustic wave penetration through the skull. These corrections allowed comparing MR-CA concentrations among animals.ResultsLow power, pulsed, ultrasound for localized disruption of BBB

  (C * ) measured at the ultrasound focal point were plotted as a function of the time elapsed between BBB disruption and MR-CA injection (Figure 4-C).

  e. closed) as a harmonic oscillator damped by fluid friction. The gap diameter as a function of the time was therefore expressed as d(t) = d 0 e -kt , where d 0 is the gap diameter just after BBB disruption and k a time constant representative of individual gap closure speed.Next, the initial size distribution of gaps was assumed to be a hemi-Gaussian function centered on 0 with a standard deviation σ 0 . These two conditions imply a distribution of gap sizes as a function of time expressed as follows:

  (n=28 rats) were used to fit this theoretical BBB closure model and determine the two parameters: a characteristic gap diameter σ 0 = 1.54 nm, and BBB closure rate k = 1.54×10 -5 s -1 .

  (d H ~ 10 nm), recombinant adeno-associated virus (rAAV) vectors for gene therapy (d H ~ 20 nm), monoclonal antibodies such as anti-β-amyloid (d H of few nm), and vectorized short interfering RNA (siRNA) used for gene expression modulation (d H ~ 50 nm).The maximum gap diameter measured in this study is probably dependent of the sonication parameters, in particular ultrasound frequency, pulse duration and duty cycle, total sonication time and size distribution of injected microbubbles (Sonovue® has a broad distribution between 1-8 µm). Here, we chose well established parameters that have been[START_REF] Choi | Noninvasive and localized blood-brain barrier disruption using focused ultrasound can be achieved at short pulse lengths and low pulse repetition frequencies[END_REF]. Several groups have studied the influence of acoustic parameters (O'Reilly et al 2011a) and microbubbles size(Samiotaki 2011) on BBB opening.

Figure 1 :

 1 Figure 1: Different steps of MRI-guided ultrasound-induced disruption of blood-brain

Figure 2 :Figure 3 :

 23 Figure 2: Ultrasound-induced localized delivery in rat brain of MRI contrast agents

Figure 4 :

 4 Figure 4: BBB closure dynamics after ultrasound induced disruption for MR-CA of

Figure 5 :

 5 Figure 5: Prediction of the time window of BBB passage for a given nanoparticle size.

Figure 3 :Figure 4 :

 34 Figure 3: Comparison of T1-weighted (T1w) and T1 mapping strategies to assess the amount of Gadolinium (Gd)-chelate crossing the BBB. Dotarem® was injected at different doses just after ultrasound-induced BBB disruption in the right thalamus, followed by the acquisition of one T1 mapping sequence (to produce a concentration map) and one T1w sequence. A. Overlay of CA concentration maps on the corresponding anatomical T2w images after injection of increasing doses of Dotarem® (375 , 750 and 1500 µmol/kg, n = 1 rat per dose). Purple arrows indicate the ultrasound focal spot where the BBB was disrupted. B. CA concentration as estimated using T1 mapping strategy at the ultrasound focal point showing linear dependence on the injected dose. C. T1w signal of the CA at the ultrasound focal point as a function of the injected dose showing a saturation effect. 124x112mm (300 x 300 DPI)

Table 1 )

 1 and verified their passage across BBB at the ultrasound focal point. The CA were intravenously

	C o n
	f i d
	e n t
	i a l :
	F
	o r
	R e v i e w
	O n l y

  quantification of the amount of MR-CA delivered to brain parenchyma in a rat model. The in designing drug delivery across the BBB, namely, the largest deliverable molecular size of the agent, and the time window for delivery of an agent of a given molecular size, were estimated for the first time.

	C o n o n C o n C o n C o n C n o two key parameters C Titles and legends to figures
	f i d f i d f i d f i d f i d f i d
	e n t e n t e n t e n t e n t e n t
	i a l : i a l : i a l : i a l : i a l : i a l :
	F F F F F F
	o r o r o r o r o r o r
	R e v i e w R e v i e w R e v i e w R e v i e w R e v i e w R e v i e w
	O n l y O n l y O n l y O n l y O n l y O n l y
	20 25

Table 1 : Principal features of the paramagnetic and superparamagnetic nanoparticle contrast agents (CA) used in this study. Hydrodynamic

 1

	1/2 (the time after

Journal of Cerebral Blood Flow and Metabolism

Acknowledgements

Authors would like to thank M. Leopoldie for providing animal care, B. Bresson for RF coil support and Dr. A. Alahari for editing this manuscript.

Disclosure/Conflict of interest

; CEA/DSV/I2BM/NeuroSpin, LRMN Larrat, Benoit; CEA/DSV/Iseult/Inumac French-German project.

Authors declare no conflict of interest.