Seongkoo Kang

Kyle G Reeves

Toshinari Koketsu

Jiwei Ma

Olaf J Borkiewicz

Peter Strasser

Alexandre Ponrouch

Damien Dambournet

Multivalent Mg 2+ , Zn 2+ and Ca 2+ Ion Intercalation Chemistry in a Disordered Layered Structure

Keywords: Disordered structure, lepidocrocite, perovskite, defects, multivalent, rechargeable batteries

The development of practical multivalent-ion batteries critically depends on the identification of suitable positive electrode materials. To gain a better understanding of the intercalation chemistry of multivalent ions, model frameworks can be used to study the distinct specificities of possible multivalent ions, thus expanding our knowledge on the emerging "Beyond Li battery" technology. Here, we compare the intercalation chemistry of Mg 2+ , Zn 2+ and Ca 2+ ions into a disordered layered-type structure featuring water interlayers and cationic vacancies as possible host sites. The thermodynamics of cation-inserted reactions performed on the model structure indicated that these reactions are thermodynamically favourable with Zn 2+ being the least stable ion. Galvanostatic measurements confirmed that the structure is inactive toward Zn 2+ intercalation while Mg 2+ can be reversibly inserted (0.37 Mg 2+ per formula unit) with minor changes of the atomic arrangement as demonstrated by pair distribution function analysis. Moreover, we demonstrate that non-solvated Mg 2+ was intercalated in the structure. Finally, the intercalation of Ca 2+ performed at 100 °C with Ca(BF 4 ) 2 in propylene carbonate, induced the collapse of the layered structure releasing water molecules that contribute to the degradation of the electrolyte as revealed by the presence of CaF 2 at the electrode level. The decomposition of the structure led to the formation of an electrochemically active phase featuring strong long-range disorder, yet short-range order close to that found in perovskite structures, particularly with corner-shared TiO 6 octahedra. We, hence, hypothesize that defective CaTiO 3based perovskite could be explored as viable cathode materials for rechargeable Ca-based batteries.

The modern-day search for new materials to be used in rechargeable batteries has focused on various components of the devices including the electrolyte and compositions of the positive and negative electrodes. An additional consideration is the ion that travels between electrodes. Given the limited amount of available lithium, previous work has explored the possibility of using other monovalent ions such as Na + or K + instead. [START_REF] Kubota | Towards K-Ion and Na-Ion Batteries as "Beyond Li-Ion[END_REF] An alternative approach takes advantage of the earth crust's abundant multivalent (MV) ions such as Mg 2+ and Ca 2+ , or Zn 2+ although it is barely more abundant than Li + , as the charge-carrying species. [START_REF] Muldoon | Quest for Nonaqueous Multivalent Secondary Batteries: Magnesium and Beyond[END_REF][START_REF] Canepa | Odyssey of Multivalent Cathode Materials: Open Questions and Future Challenges[END_REF][START_REF] Aurbach | Prototype Systems for Rechargeable Magnesium Batteries[END_REF][START_REF] Ponrouch | Towards a Calcium-Based Rechargeable Battery[END_REF][START_REF] Arroyo-De Dompablo | Challenges, and Prospects of Calcium Batteries[END_REF][START_REF] Konarov | Present and Future Perspective on Electrode Materials for Rechargeable Zinc-Ion Batteries[END_REF][START_REF] Han | Origin of Electrochemical, Structural, and Transport Properties in Nonaqueous Zinc Electrolytes[END_REF] In a system which uses a metal anode, multivalent ions take advantage of a multi-electron processes leading to a corresponding increase in the volumetric energy density (Table 1). The proper plating/stripping of the metal in a selected electrolyte without the formation of dendrite is a mandatory step toward the realization of beyond lithium-ion batteries with enhanced safety. However, while the energy density of a battery can, in principle, be significantly enhanced using Mg, Ca or Zn metal anodes, the development of positive electrode materials that can accommodate these ions is still in its early stages. [START_REF] Rong | Materials Design Rules for Multivalent Ion Mobility in Intercalation Structures[END_REF][START_REF] Canepa | Odyssey of Multivalent Cathode Materials: Open Questions and Future Challenges[END_REF] The difficulty lies in the slow solid-state diffusion of MV ions within the structure caused by the strong Coulombic interactions with the anionic sublattice. To account for this, the polarizing power (ratio between the charge and the square of the ionic radius: Z/r 2 ) reflects the ability of a cation to attract and distort the electron cloud of a nearby anion (Table 1). Considering a six-fold coordination number (CN), the ranking goes Li + < Ca 2+ < Zn 2+ < Mg 2+ indicating that the greater challenges may be expected in developing beyond-lithium ion rechargeable batteries. This trend, however, does not take into account the specific behaviour of these ions within electrode structures, an environment unlike bulk solvation in which ions are confined. These potential new generations of batteries require more investigation to better understand the stability of each ion species when intercalated into the electrode material. [START_REF] Rong | Materials Design Rules for Multivalent Ion Mobility in Intercalation Structures[END_REF] Model frameworks such as layered structures are particularly interesting to explore due to the ability of these materials to accommodate ions with a range of sizes and the possible variability in the coordination environment for intercalated species. [START_REF] Cui | Recent Advances in the Rational Design and Synthesis of Two-Dimensional Materials for Multivalent Ion Batteries[END_REF] Moreover, highly defective structures were shown to favour the diffusion of MV ions. [START_REF] Koketsu | Reversible Magnesium and Aluminium Ions Insertion in Cation-Deficient Anatase TiO2[END_REF][START_REF] Orikasa | Noncrystalline Nanocomposites as a Remedy for the Low Diffusivity of Multivalent Ions in Battery Cathodes[END_REF] Herein, we embark on comparing the intercalation chemistry of MV ions within a disordered layered structure featuring interlayer water molecules and cationic vacancies both potentially acting as host sites. DFT-calculations on simplified structural model were primarily used to access the thermodynamic properties and to provide structural insights into the properties of inserted cations highlighting noticeable differences. Although the comparison between ions is constrained by the specific electrochemical features of each cell (e.g. the electrolyte associated with each multivalent ion), this work provides fundamental insights into different intercalation chemistries.

EXPERIMENTAL SECTION

Synthesis. X-ray amorphous sample with a local lepidocrocite-type layered structure containing a high concentration of titanium vacancies was synthesized by a sol-gel process as previously reported. [START_REF] Ma | Layered Lepidocrocite Type Structure Isolated by Revisiting the Sol-Gel Chemistry of Anatase TiO2: A New Anode Material for Batteries[END_REF] Briefly, a solution containing 4 mL of titanium (IV) isopropoxide (> 97 %, Sigma-Aldrich), 25.19 mL of isopropanol (Alfa Aesar) and 0.81 mL of ultrapure water, was placed in a Teflon-lined container sealed in the autoclave. The mixture underwent thermally activated precipitation at 90 °C for 12 h. The obtained white precipitate was washed several times with ethanol (> 96 %), collected by centrifugation and further dried at 80 °C overnight. The determination of the chemical formula led to Ti 1.5  0.5 O 2 (OH) 2 •nH 2 O ( = Ti 4+ vacancy, n~0.55) where 25% Ti 4+ replaced by vacancies to compensate the charge difference due to the substitution of divalent O 2-by monovalent OH -groups.

Electrochemistry.

The electrochemical activity of Ti 1.5  0.5 O 2 (OH) 2 •nH 2 O versus Mg 2+ , Ca 2+ and Zn 2+ was studied with galvanostatic charge-discharge experiments. The electrodes were prepared by mixing different proportions of active materials, conductive carbon (Super P, Timcal) and polyvinylidene fluoride (PVDF) as the binder. Galvanostatic cycling measurements were done in three-electrode Swagelok cells using corresponding metal wire and foil as reference and counter electrode and suitable electrolyte to electrochemically strip and deposit the metal. Experiments were performed at room temperature for Mg and Zn, and at 100 °C for Ca cells. The current density was 10 mA.g -1 for Mg and Zn and 5 mA.g - 1 for Ca. Detailed experimental conditions can be found in Table 2. Ex situ Synchrotron X-ray scattering. Electrochemically reduced/oxidized electrodes were washed with tetrahydrofuran or dimethyl carbonate and scrapped from the current collector and packed in kapton capillaries. Operations were performed in an Ar filled glove box. High energy X-ray (λ = 0.2128 Å) synchrotron diffraction data were collected at the 11-ID-B beamline at the Advanced Photon Source at Argonne National Laboratory. 1D diffraction data were obtained by integration of 2D total scattering data with Fit2D. [START_REF] Hammersley | Two-Dimensional Detector Software: From Real Detector to Idealised Image or Two-Theta Scan[END_REF] Pair Distribution Functions (PDF) were extracted using PDFgetX3. [START_REF] Juhás | PDFgetX3 : A Rapid and Highly Automatable Program for Processing Powder Diffraction Data into Total Scattering Pair Distribution Functions[END_REF] The PDF refinements were performed using PDFgui [START_REF] Farrow | PDFfit2 and PDFgui: Computer Programs for Studying Nanostructure in Crystals[END_REF] . Energy dispersive X-ray spectroscopy. EDX analysis was performed using a JEOL 2100 Plus transmission electron microscopy LaB6 200 kV. Density function theory calculations. A lamellar lepidocrocite structure was simulated using a supercell made of two titanate layers (4x4x1 supercell) and two layers of interlayer water with an overall composition Ti 63  1 O 124 (OH) 4 •32H 2 O where the  represents a single Ti 4+ vacancy whose charge is compensated by the four hydroxides. The spacing between titanate layers was fixed at 11.5 Å to be consistent with existing literature [START_REF] Yuan | The Swelling Transition of Lepidocrocite-Type Protonated Layered Titanates into Anatase Under Hydrothermal Treatment[END_REF] . The simulation cell was computed using periodic boundary conditions. The choice to simplify the system was made to control for the role of the vacancy and to eliminate significant contributions from disorder in the structure. From this approach, this work will investigate the behaviour of intercalated multivalent ions (e.g. Mg 2+ , Ca 2+ and Zn 2+ ) and the structure and thermodynamics of stable structures in the vacancy and inter-layer positions. Density functional theory (DFT) calculations were performed using the CP2K code via the Quickstep algorithm. [START_REF] Hutter | Atomistic Simulations of Condensed Matter Systems[END_REF][START_REF] Vandevondele | Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach[END_REF] Electrons were described using the PBE exchange-correlation function with a planewave cut-off of 400 Rydberg. DZVP-MOLOPT-SR-GTH basis sets were used along with Goedecker-Teter-Hutter pseudopotentials [START_REF] Vandevondele | Gaussian Basis Sets for Accurate Calculations on Molecular Systems in Gas and Condensed Phases[END_REF][START_REF] Goedecker | Separable Dual-Space Gaussian Pseudopotentials[END_REF] . Preliminary tests using the Hubbard-U correction (U = 2.5 eV) were performed to understand the influence of correlation effects on energies of formations. A table summarizing these results can be found in the SI. DFT calculations in the article did not used the Hubbard-U correction as the differences in the computed energies of formation were on the order of 10meV/formula unit. Given the size of the simulation cell, only the Γ-point was sampled. Systems are computed with a total neutral charge to simulate the state of the electrode following a discharge process. Energy of Formation. Based on the Wyckoff positions, the energy of formation was computed for two different sites within the lepidocrocite structure: the ion in the water interlayer and the ion in the Ti 4+ vacancy as shown in Figure 1B and C. To determine the energy of formation of the ion within the lepidocrocite structure, we relax the structure for the ion in each of the two previously mentioned positions allowing for both the ion and the electrode atoms to freely move. When this relaxed structure is determined, the total energy of the system is then recomputed for this arrangement of atoms but without the intercalated ion. The energy of formation is then computed via the formula:

E E i E E i
where the first term represents the combined system and the terms in parentheses represent the sum of energies of the ion, and titanate systems computed independently.

RESULTS AND DISCUSSION

Thermodynamic and analysis of the coordination environment of inserted-cations

The thermodynamic and coordination environment of each cation was probed by using DFT-calculations performed on a simplified structure, i.e., a supercell Ti 63

 1 O 124 (OH) 4 •32H 2 O
where the  represents a single Ti 4+ vacancy whose charge is compensated by the four hydroxides, as shown in Figure 1.

In both cases of the interlayer Wyckoff site (Figure 1B) and the cationic vacancy (Figure 1C), formation energies indicate that the intercalation reaction are thermodynamically favourable (Figure 2). We, however, note that zinc is the least stable ion within the two sites as compare to calcium and magnesium. We detailed below the structural features for each cation. Magnesium. After completing the geometry relaxation, the magnesium ion is four-fold coordinated with water molecules (Figure 3A). The ion remains in the centre of the interlayer suggesting that the structure can remain stable upon Mg 2+ insertion. Within the vacancy (Figure 3D), Mg 2+ ion is 6-fold coordinated with oxide anions with an average bond length of 2.06 Å. The ion is centred within the vacancy in almost exactly in the position one would expect for Ti 4+ . As a result, there is only a small amount of distortion associated with the insertion of the ion into the vacancy site. Calcium. The relaxed structure computed for the interlayer calcium ion shows that the ion moved away from the centre of the interlayer. As a result, it is solvated by three water molecules and coordinated with a bridging oxygen from the titanate layer (Figure 3B). Such a configuration suggests that calcium insertion might destabilize the structure by moving water molecules away for the interlayer space. Within the vacancy, the ion forms four shorter bonds (~2.22 Å) and two longer bonds (~2.48 Å) (Figure 3E) effective coordination number 24 of 4.4. This is due to the position of the calcium ion with respect to the centre of the vacancy. Unlike the Mg 2+ , which falls towards the centre of the vacancy, the larger size of Ca 2+ induced a position farther away from the vacancy in the direction of the interlayer such that it is coordinated principally by the four outermost oxygen atoms. This destabilization as result of the distortion to the titanate structure has previously been observed in the case of K + in the lepidocrocite structure. [START_REF] Reeves | Insights into Li+, Na+, and K+ Intercalation in Lepidocrocite-Type Layered TiO2 Structures[END_REF] Zinc intercalation. When intercalated into the interlayer, Zn 2+ ion is three-fold coordinated by water molecules (Figure 3C).

Similar to Mg 2+ ion, zinc remains at the centre of the interlayer with no interactions with the titanate layer. When the ion is intercalated into the vacancy site, it sits higher than the centre of the void leading to a four-fold coordinated (Figure 3F), which is its preferred coordination number [START_REF] Brown | What Factors Determine Cation Coordination Numbers?[END_REF] . The four bonds formed with the Zn 2+ have an average length of 2.07 Å with the two remaining oxygen atoms falling at 2.82/2.73 Å away from the ion. Such a configuration suggests an instability of Zn 2+ to sit fully within the vacancy.

Electrochemistry. The electrode material that was used here consists of an amorphous compound featuring a local structure previously described using the layered lepidocrocite network.

Disorder was suggested to be due to the large content of vacancies stabilized by the partial substitution of divalent O 2-by monovalent OH -groups, i.e., Ti 1.5  0.5 O 2 (OH) 2 •nH 2 O ( = Ti 4+ vacancy, n~0.55). [START_REF] Ma | Layered Lepidocrocite Type Structure Isolated by Revisiting the Sol-Gel Chemistry of Anatase TiO2: A New Anode Material for Batteries[END_REF] The electrochemical properties of the lepidocrocite-type structure with respect to MV ions were assessed by employing the corresponding metal as the negative electrode. Moreover, specific electrolytes were chosen for each system (Table 2). Figure 4 gathers galvanostatic charge/discharge curves obtained for the first three cycles. In the case of Mg 2+ , the all-phenyl complex (APC) Mg 2 Cl 3 + AlPh 2 Cl 2 -dissolved in tetrahydrofuran (THF) was used as the electrolyte. The voltage-composition (Figure 4a) curves measured at RT under 10 mA.g -1 , features a smooth profile suggesting a solid-solution behaviour. The reversible capacity was 135 mAh g -1 which corresponds to the intercalation of 0.37 Mg 2+ per formula unit. Energy-dispersive x-ray spectroscopy was performed on a discharged electrode (Figure S1a) confirmed the insertion of naked Mg 2+ with negligible amount of chloride (Cl/Ti = 0.02) ruling out the insertion of solvated species such as magnesium monochloride cations [START_REF] Yoo | Fast Kinetics of Magnesium Monochloride Cations in Interlayer-Expanded Titanium Disulfide for Magnesium Rechargeable Batteries[END_REF] . In the case of Zn 2+ ions, the structure was shown to be electrochemically inactive (Figure 4b). The inability of the network to accommodate Zn 2+ was confirmed by testing other type of solvent such as dimethyl sulfoxide. In the latter case, we measured very low capacities. The inability of Zn 2+ to be inserted in the layered structure agrees with the thermodynamic reinforcing the observation of Zn as the least stable cation in the structure. For Ca 2+ intercalation, the selected electrolyte which was 0.45 M Ca(BF 4 ) 2 / EC:PC (50:50 wt%) requires to operate at 100 °C to provide an efficient plating/stripping of calcium. [START_REF] Ponrouch | Towards a Calcium-Based Rechargeable Battery[END_REF] Note that the thermal analysis showed that the structure is stable at this temperature. [START_REF] Ma | Layered Lepidocrocite Type Structure Isolated by Revisiting the Sol-Gel Chemistry of Anatase TiO2: A New Anode Material for Batteries[END_REF] The first discharge led to about 290 mAh.g -1 (Figure 4c). EDX analysis performed on a discharge electrode indicated the presence of fluorine pointing to the decomposition of the boron fluoride groups from Ca(BF 4 ) 2 (Figure S1). The subsequent charge capacity reached 134 mAh.g -1 , corresponding to about 0.35 Ca 2+ reversibly inserted per formula with an operating voltage of about 2.25V vs. Ca pseudo reference electrode. Large irreversible capacity is observed during the first few cycles assigned to the decomposition of the electrolyte and formation of a passivation layer containing calcium fluoride. Coulombic efficiency above 98% is reached only after 5 cycles (Figure S2), which is probably due to the progressive decomposition of the electrolyte and formation of CaF 2 -based layer. Intercalation mechanisms. The intercalation mechanisms of Mg 2+ and Ca 2+ within the lepidocrocite structure were investigated using the pair distribution function. Such a method obtained by the Fourier transformation of high-energy x-ray data, enables to get structural information of x-ray amorphous compounds. [START_REF] Billinge | Beyond Crystallography: The Study of Disorder, Nanocrystallinity and Crystallographically Challenged Materials with Pair Distribution Functions[END_REF] Figure 5 gathered the PDFs of the pristine, discharged and charged electrodes vs. Mg 2+ /Mg. The PDF of the pristine electrode presents peaks whose intensities gradually decrease and vanish at inter-atomic distances that are greater than 10 Å, which is a characteristic of x-ray amorphous compounds. The insertion and subsequent de-insertion of Mg 2+ produced minor changes to the PDFs suggesting that the local structure is maintained upon (de)-intercalation. Upon Mg 2+ insertion, the intensity of the first peak characteristic of Ti-O(H) bond distances decreased while a small shoulder at 2.16 Å appears, which is assigned to reduced Ti 3+ -O(H) bonds. This implies that the redox activity at titanium sites is a consequence of the Ti 4+ /T i3+ redox couple. Moreover, the second peak at 3.1 Å can be assigned to edge-shared TiTi bonds and the third peak at 3.65 Å can be attributed to TiO distances broaden upon Mg 2+ insertion, which can be due to the appearance of reduced Ti 3+ -O(H) and Mg-O(H) bonds. Small features can be observed at higher inter-atomic distances, i.e., ~9 Å, suggesting that the Mg 2+ intercalation produced changes up to the intermediate range-order, ruling out capacitive contributions. We further attempted to refine the PDF data using a disordered lepidocrocite structure (Figure S3). The results, however, was not conclusive, particularly in discriminating any sort of site preferences for Mg-ions. It can be postulated that Mg 2+ can be inserted at both sites, and the populations in each site is a result of the relative stability the ion experiences either via the solvation of water molecules or within the titanate layer. We also note that Mg 2+ insertion take place without solvated chloride-based species. Upon subsequent Mg 2+ de-intercalation, the peaks match those of the pristine showing that Mg 2+ insertion proceeds in a reversible manner. We note that even at the intermediate rangeorder, the peaks match their initial position. The network stability upon Mg 2+ intercalation/de-intercalation primary suggests that water molecules remain in the interlayer space and does not shuttle into the electrolyte. The precise role of interlayer water molecules during Mg 2+ intercalation and deintercalation in 2D hydrated structure is still challenging to establish. [START_REF] Tepavcevic | Nanostructured Layered Cathode for Rechargeable Mg-Ion Batteries[END_REF][START_REF] Sa | Structural Evolution of Reversible Mg Insertion into a Bilayer Structure of V2O5•nH2O Xerogel Material[END_REF][START_REF] Ji | Water-Activated VOPO4 for Magnesium Ion Batteries[END_REF][START_REF] Sai Gautam | Role of Structural H2O in Intercalation Electrodes: The Case of Mg in Nanocrystalline Xerogel-V2O5[END_REF][START_REF] Tan | Crystal Regulation towards Rechargeable Magnesium Battery Cathode Materials[END_REF] In particular, the structure stability upon extended cycling would need a deeper investigation. To better understand Mg 2+ intercalation within the lepidocrocite structure, we compared the energy of formation for several positions of the ion throughout the structure of the computed supercell Ti 63  1 O 124 (OH) 4 •32H 2 O. As shown by a dotted line in Figure 6, a path perpendicular to the titanate layer passing between water interlayers via the titanium vacancy was considered. At each position of the ion along this path, the structure was allowed to relax to find a local energy minimum. Thus, the profile of energies in Figure 6 should not be taken as the diffusion pathway of ions, but a series of local minima computed to provide insight into the interaction between the ions and the host material. We compare the two profiles for Li + and Mg 2+ , highlighting three positions along the ion path: in the titanium vacancy (Position A), just outside the vacancy (Position B), and in the interlayer (Position C). To better compare the energy of formation profiles between ions, Figure 6 shows the difference in the energy of formation (ΔE) with respect to the lowest computer energy formation along the path. For each ion, the lowest energy of formation was determined to be in the interlayer position. It should be noted that an increase in the profile represents relative stability amongst the varying positions. Energies of formation remain favourable, as was shown in Figure 2. Moreover, we expected the diffusion of Mg 2+ to be more favourable in the disordered material, particularly when approaching the titanate layer, i.e., when the formation energy increases drastically. The profiles at the three noted positions, do not vary significantly, although most notable at position B (~0.72eV). This position, corresponding to the ion just outside of the vacancy, looks structurally similar between the two ions. In each case, the ion prefers to move out of the plane between the two hydroxyl groups, minimizing the amount of distortion that is experienced within the structure. The distance between the hydroxyl groups in the relaxed structures were 3.69Å and 3.67 Å for Li + and Mg 2+ respectively. The similarity in these distances suggests that the difference in energy is not likely due to differing amounts of distortion induced by the ions themselves, and we thus conclude that the polarizing power of the ion is likely the main driving force in the energetic penalty of the magnesium ion. As opposed to Mg 2+ , the insertion of Ca 2+ induced a drastic structural change. X-ray diffraction analysis (Figure S1) of the discharged electrode revealed the presence of a crystallized compound attributed to CaF 2 (S.G: Fm-3m) and weak peaks tentatively assigned to CaB 3 O 4 (OH) 4 H 2 O 34 . Both phases indicate the degradation of the electrolyte. The PDF of the discharged electrode (Figure 7a) showed peaks beyond 10 Å, which can be attributed to the presence of newly formed crystallized phases. Moreover, for low r region, the peak located at 3.1 Å assigned to edge-shared TiTi bond distances disappeared, suggesting a phase change with new atoms connectivity's formed upon Ca 2+ insertion. The disappearance of edge-shared TiTi connectivities indicated that the layer structure was not maintained and therefore, that water molecules were released in the electrolyte. This is supported by the large content of CaF 2 observed by XRD, whose formation can be promoted by the hydrolysis of BF 4 -groups releasing fluoride ions [START_REF] Freire | Hydrolysis of Tetrafluoroborate and Hexafluorophosphate Counter Ions in Imidazolium-Based Ionic Liquids[END_REF] . To understand the structural changes induce by Ca 2+ insertion, we attempted to refine the PDF data with the following procedure: (i) we refined the data by excluding the short-range order (i.e., for r > 7 Å). For 7-30 Å, the PDF was well reproduced by using CaF 2 as a structural model. Subsequently, we fixed all the parameters and refined the data from 1 to 20 Å (Figure 7b), which revealed inter-atomic distances that can be assigned to the newly formed phase. Such a phase features only a short-range order. We tentatively assigned the most intense peaks at 1.94 Å to Ti-O bonds, at 2.3 Å to Ca-O bonds and at 3.8 Å to corner-shared TiTi. Considering the absence of edge-shared TiTi connectivities, we deduced from the ICSD database that the newly formed phase has a local perovskitelike structure as found in CaTiO 3 [START_REF] Sasaki | Orthorhombic Perovskite CaTiO3 and CdTiO3: Structure and Space Group[END_REF] , a structure where Ca 2+ is 7-fold coordinated (Inset Figure 7d), which is close to its preferred 8-fold coordination number [START_REF] Rong | Materials Design Rules for Multivalent Ion Mobility in Intercalation Structures[END_REF][START_REF] Brown | What Factors Determine Cation Coordination Numbers?[END_REF] . The calculated PDF of CaTiO 3 (Figure S2) is, however, characterized by an intense peak at 3.2 Å assigned to CaTi connectivities. The weak presence of this peak in the newly formed phase suggested an under stoichiometry in Ca 2+ . The exact chemical formula is, however, not trivial to determine. We simulated the PDF of CaTiO 3 (Figure 7c) featuring low content of Ca (Ca occupancy = 0.34) and found that it matches the experimental data with a PDF scattering domain length of around 8 Å highlighting strong disorder. The newly formed phase enabled to insert Ca 2+ in a reversible manner. Our work confirms that disordered perovskite-based structure could be seen as an attractive framework to enable practical Ca-batteries. 37

CONCLUSION

We investigated the electrochemical properties vs. multivalent ions (Mg 2+ , Zn 2+ , Ca 2+ ) of a disordered layered type structure featuring water-containing interlayer and titanium vacancies (Ti 1.5  0.5 O 2 (OH) 2 •nH 2 O). First-principle calculations were primary performed on a simplified supercell Ti 63  1 O 124 (OH) 4 •32H 2 O to assess the thermodynamics of ions intercalation reaction within either the interlayer space or titanium vacancy. Despite favourable formation energies, experimental data showed drastic differences between ions due to their intrinsic properties (size, polarizing power) and operating conditions (temperature, electrolyte). Galvanostatic experiments showed that Zn 2+ cannot be intercalated within the lepidocrocite-type structure possibly due to unfavorable arrangement within the vacant sites as revealed by firstprinciple calculations. In the case of magnesium, the structure is able to reversibly intercalated 0.37 Mg 2+ per chemical formula. Pair distribution function computed for the discharged and charged electrodes demonstrated that the layered framework is maintained. We suggest that this is due to only slight distortions to the titanate structure as predicted by geometry relaxations and due to interlayer water maintaining its role of structural stabilization during the insertion of Mg 2+ . Moreover, unsolvated Mg 2+ ions are inserted as demonstrated by EDX analysis showing the absence of chlorine, ruling out the insertion of magnesium chlorinated species. The electrochemical insertion of Ca 2+ induced the collapse of the layered structure with the release of the interlayer water molecules that contributed to the degradation of the electrolyte as shown by the formation of CaF 2 . The refinement of the pair distribution function measured on the discharged electrode indicated that the newly formed phase, which is electrochemically active, featured strong disorder and short-range order close to those found in perovskite structure, particularly with corner-shared TiO 6 octahedra. Perspectives in synthetizing defective CaTiO 3 -based perovskite could be explored to identify viable cathode materials for rechargeable Ca-based batteries.

ASSOCIATED CONTENT Supporting Information

Figure 1 .

 1 Figure 1. (A) The model lepidocrocite system considered here contains two titanate and two interlayer water layers. A single vacancy is introduced into one of the two titanate layers by removing a Ti atom and balancing the charge with four protons. Multivalent ions were inserted into one of two Wyckoff sites: the interlayer or the vacancy position (B) and (C) respectively, where the position of the multivalent ions is shown in green.

Figure 2 .

 2 Figure 2. Energies of formation (meV/formula unit) for multivalent ions located in vacancy and interlayer Wyckoff positions in lepidocrocite-type TiO 2 .

Figure 3 .

 3 Figure 3. The relaxed geometries for Mg, Ca, and Zn ions in the water interlayer (A-C) and in the Ti 4+ vacancy (D-F). Bond lengths are shown in units of Angstroms.

Figure 4 .

 4 Figure 4. Galvanostatic discharge-charge curves in a) Mg, b) Zn and c) Ca cells. The used current densities were 10 mA.g -1 for Mg 2+ and Zn 2+ , and 5 mA.g -1 for Ca 2+ .

Figure 5 .

 5 Figure 5. Ex-situ X-ray atomic pair distribution function of pristine, fully discharged and fully charged sample vs. Mg 2+ /Mg, a) extended and b) zoomed view.

Figure 6 .

 6 Figure 6. Change in the energy of formation for Li + (green) and Mg 2+ (orange) ions in the lepidocrocite structure. Relaxed geometries are visualized at the vacancy structure (A), having just left the titanium vacancy (B), and the water interlayer (C). The energy of formation profile corresponding to a path between interlayers, shown as a dotted black line.

Figure 7 .

 7 Figure 7. a) X-ray PDFs of pristine and fully discharged electrode vs. Ca 2+ /Ca. b) refinement of the PDF data using CaF 2 (S.G. Fm-3m) as a structural model and c) using CaF 2 and Ca x TiO 3 (S.G. Pbnm) perovskite structure, d) Polyhedral representation of the perovskite structure.

Table 1 .

 1 Characteristics of different charge carriers.

	Characteristics		Li +	Ca 2+	Zn 2+	Mg 2+
	Earth crust abundance	20	41,000 75	23,300
	(ppm)
	Volumetric capacity	2,061	2,073	5,854	3,832
	(mAh/cm 3 )
	Shannon's ionic radii	0.76	1.00	0.74	0.72
	(Å) (CN*=VI)				0.60**
	Polarizing	power	1.73	2.00	3.65	3.85
	(Z/r 2 ) (Å -2 )				5.55**

*CN: Coordination number ** CN of Zn 2+ in ZnO is IV.

Table 2 .

 2 Experimental conditions for electrochemical measurements.

	Counter	Electrode formulation:	Electrolyte 5,8,15
	electrode	AM, C, PVDF
		Weight %
	Mg	80 wt%, 10 wt%, 10	0.2 M 2PhMgCl-
		wt%	AlCl 3 / THF
	Zn	70 wt%, 20 wt%, 10	0.5 M Zn(TFSI) 2 /
		wt%	acetonitrile
	Ca	80 wt%, 10 wt%, 10	0.45 M Ca(BF 4 ) 2 /
		wt%	EC:PC	(50:50
			wt%)

ACKNOWLEDGMENT

The work done at the Advanced Photon Source, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science by Argonne National Laboratory, was supported by the U.S. DOE under contract no. DE-AC02-06CH11357. A. Ponrouch is grateful to the Spanish Ministry for Economy, Industry and Competitiveness Severo Ochoa Programme for Centres of Excellence in R&D (SEV-2015-0496). We thank IMPC FR2482 for SEM-FEG and TEM instrumentation funded by CNRS, Sorbonne Université and C'Nano projets of Region Ile-de-France.

AUTHOR INFORMATION

Corresponding Author * damien.dambournet@sorbonne-universite.fr Author Contributions ‡ hese authors contributed equally.

Funding Sources

The research leading to these results has received funding from the French National Research Agency under Idex@Sorbonne University for the Future Investments program (no. ANR-11-IDEX-0004-02). Funding from the European Union's Horizon 2020 research is acknowledged: European Research Council (ERC-2016-STG, CAMBAT grant agreement No 715087).


Nazim Haouchine

Parikshit Juvekar

William M Wells

Stephane Cotin

Alexandra Golby

Sarah Frisken

Deformation Aware Augmented Reality for Craniotomy using 3D/2D Non-rigid Registration of Cortical Vessels

Keywords: 3D/2D Non-rigid registration, Physics-based Modelling, Augmented Reality, Image-guided Neurosurgery, Shape-from-template

Intra-operative brain shift is a well-known phenomenon that describes non-rigid deformation of brain tissues due to gravity and loss of cerebrospinal fluid among other phenomena. This has a negative influence on surgical outcome that is often based on pre-operative planning where the brain shift is not considered. We present a novel brain-shift aware Augmented Reality method to align pre-operative 3D data onto the deformed brain surface viewed through a surgical microscope. We formulate our non-rigid registration as a Shape-from-Template problem. A pre-operative 3D wire-like deformable model is registered onto a single 2D image of the cortical vessels, which is automatically segmented. This 3D/2D registration drives the underlying brain structures, such as tumors, and compensates for the brain shift in sub-cortical regions. We evaluated our approach on simulated and real data composed of 6 patients. It achieved good quantitative and qualitative results making it suitable for neurosurgical guidance.

Introduction

Brain shift is a well-known intra-operative phenomenon that occurs during neurosurgical procedures. This phenomenon consists of deformation of the brain that changes the location of structures of interest from their locations in preoperative imaging [START_REF] Bayer | Intraoperative imaging modalities and compensation for brain shift in tumor resection surgery[END_REF] [START_REF] Miga | Clinical evaluation of a model-updated image-guidance approach to brain shift compensation: experience in 16 cases[END_REF]. Neurosurgical procedures are often based on preoperative planning. Therefore, estimating intra-operative brain shift is important since it may considerably impact the surgical outcome. Many approaches have been investigated to compensate for brain shift, either using additional intra-operative imaging data [START_REF] Kuhnt | Brain shift compensation and neurosurgical image fusion using intraoperative mri: Current status and future challenges[END_REF] [20] [START_REF] Ji | Mutual-informationbased image to patient re-registration using intraoperative ultrasound in imageguided neurosurgery[END_REF] [START_REF] Rivaz | Deformable registration of preoperative mr, pre-resection ultrasound, and post-resection ultrasound images of neurosurgery[END_REF] or advanced brain models to predict intra-operative outcome pre-operatively [START_REF] Bilger | Biomechanical simulation of electrode migration for deep brain stimulation[END_REF] [8] [START_REF] Essert | Automatic computation of electrode trajectories for deep brain stimulation: A hybrid symbolic and numerical approach[END_REF]. The latter approach has the advantage of being hardware-independent and does not involve additional imaging in the operating room. However, even sophisticated brain models have difficulty modeling events that occur during the procedure. Additional data acquisition is then necessary to obtain a precise brain shift estimation. Various types of imaging techniques have been proposed to acquire intra-operative information such as ultra-sound [START_REF] Ji | Mutual-informationbased image to patient re-registration using intraoperative ultrasound in imageguided neurosurgery[END_REF] [21] [START_REF] Reinertsen | Intraoperative correction of brain-shift[END_REF], Cone-Beam Computed Tomography [START_REF] Pereira | Volumetric Measurements of Brain Shift Using Intraoperative Cone-Beam Computed Tomography: Preliminary Study[END_REF] or intra-operative MRI [START_REF] Kuhnt | Brain shift compensation and neurosurgical image fusion using intraoperative mri: Current status and future challenges[END_REF].

During a craniotomy, the cortical brain surface is revealed and can be used as an additional source of information. Filipe et al. [START_REF] Marreiros | Non-rigid deformation pipeline for compensation of superficial brain shift[END_REF] used 3 Near-Infrared cameras to capture brain surface displacement which is then registered to MRI scans using the coherent point drift method by considering vessel centerlines as strong matching features. A new deformed MRI volume is generated using a thin plate spline model. The approach presented by Luo et al. [START_REF] Luo | A comprehensive model-assisted brain shift correction approach in image-guided neurosurgery: a case study in brain swelling and subsequent sag after craniotomy[END_REF] uses an optically tracked stylus to identify cortical surface vessel features; a model-based workflow is then used to estimate brain shift correction from these features after a dural opening. Jiang et al. [START_REF] Jiang | Marker-less tracking of brain surface deformations by non-rigid registration integrating surface and vessel/sulci features[END_REF] use phase-shifted 3D measurement to capture 3D brain surface deformations. This method highlights the importance of using cortical vessels and sulci to obtain robust results. The presence of a stereo-miscroscope in the operating room makes it very convenient to deploy such methods clinically. Sun et al. [START_REF] Sun | Near realtime computer assisted surgery for brain shift correction using biomechanical models[END_REF] proposed to pre-compute brain deformations to build an atlas from a sparse set if image-points extracted from the cortical brain surface. These image-points are extracted using an optically tracked portable laser range scanner. Haouchine et al. [START_REF] Haouchine | Alignment of cortical vessels viewed through the surgical microscope with preoperative imaging to compensate for brain shift[END_REF] proposed to use a finite element model of the brain shift to propagate cortical deformation captured from the stereoscope to substructures. In a similar way Mohammadi et al. [START_REF] Mohammadi | Estimation of intraoperative brain shift by combination of stereovision and doppler ultrasound: phantom and animal model study[END_REF] proposed a projection-based stereovision process to map brain surface with a pre-operative finite element model. A pre-defined pattern is used to recover the 3D brain surface. In order to build a dense brain surface and gather more precise positions, Ji et al. [START_REF] Ji | Cortical surface shift estimation using stereovision and optical flow motion tracking via projection image registration[END_REF] proposed a stereo-based optical flow shape reconstruction. The 3D shapes are recovered at different surgical stages to obtain undeformed and deformed surfaces. These surfaces can be registered to determine the deformation of the exposed brain surface during surgery.

Contribution:

We propose a novel method to register pre-operative scans onto intra-operative images of the brain surface during the craniotomy. As shown in the pipeline of figure 1, our method uses images from a surgical microscope (or possibly a ceiling mounted camera), rather than intra-operative ultrasound or MRI which requires significant cost and time. Unlike previous methods, we rely solely on a single image to avoid tedious calibration of the stereo camera, laser range finder or optical stylus. Our method considers cortical vessels centerlines as strong features to drive the non-rigid registration. The intra-operative centerlines are automatically extracted from the image using convolutional deep neural networks, while the pre-operative centerlines are modelled as a network of linked beams capable of handling non-linear deformations. Our approach is formulated as a force-based shape-from-template problem to register the undeformed pre-operative 3D centerlines onto the deformed intra-operative 2D centerlines. This problem is solved by satisfying a combination of physical and projective constraints. We present our results through the microscope occulars using Augmented Reality view by overlaying the tumor model in the miscropcic view after accounting for the estimated brain shift deformation. 2 Method

Extracting 2D Cortical Brain Vessels

To extract the vessels from microscopic images we rely on a deep convolutional network following a typical U-Net architecture [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. The input is a single RGB image of size 256 × 256 × 3 and the final layer is a pixel-wise Softmax classifier that predicts a class label {0, 1} for each pixel. The output image is represented as a binary image of size 256 × 256 as illustrated in Figure 2. the segmentation can suffer from class imbalance because microsopic images of the cortical surface are composed of veins, arteries and parenchyma of different sizes. We thus rely on a weighted cross-entropy loss L to account for imbalanced depths [START_REF] Eigen | Predicting depth, surface normals and semantic labels with a common multi-scale convolutional architecture[END_REF], formulated as follows:

L = - 1 n i w ci •c * i •log(c i ) with c i = e zi c e zi,ci and w ci = med(f ci |c i ∈ C) f ci ( 1
)
where n is the number of pixels, C is the set of classes, c * i is the ground truth label of pixel i, c i is the probability of classifying pixel i with class c i . (z i , c i ) is the output of the response map and w ci is the class-balanced weight for class c i , with f ci being the frequencies of class c i and med(•) the median. Once segmented, the vessel centerlines are extracted from the binary image by skeletonization. We denote v, the vector of 2D points corresponding to the intra-operative centerlines.

3D Vascular Network Modelling

Let us denote u the vector of 3D vertices representing the vessel centerlines derived from MRI scans. This set of vertices is used to build a non-rigid wire-like model that represents the behaviour of cortical vessels. More precisely, the vessels are modeled with serially linked beam elements that can handle geometric non-linearities while maintaining real-time computation. This model has previously been used to simulate guide wires and catheters [START_REF] Cotin | New approaches to catheter navigation for interventional radiology simulation[END_REF]. Each beam element is delimited by two nodes each having 3 spatial and angular degrees of freedom. These node positions relate to the forces applied to them thanks to a stiffness matrix K e and account for rotational degrees of freedom through a rotational matrix R e . At each node i, the internal forces f i generated by the deformation of the structures is formulated as:

f i = i e=i-1 R e (u)K e (R e (u) T (u -u e ) -u rest ) ( 2
)
where e is the index of the two beams connected to this i th node. u i-1 , u i and u i+1 are the degrees of freedom vectors of the three nodes (respectively i -1, i, i + 1) and belong to the two beams in the global frame. u ej denotes the middle frame of the j th beam that is computed as an intermediate between the two nodes of the beam, where u rest corresponds to the degrees of freedom at rest in the local frame. The global force f emanating from the 3D vessels can be computed as f = nu j f j , with n u being the number of nodes of the vascular tree.

Force-based Shape-from-Template

A shape-from-template problem is defined as finding a 3D deformed shape from a single projected view knowing the 3D shape in rest configuration [START_REF] Bartoli | Shape-fromtemplate[END_REF]. In our case, the 3D shape at rest consists of the pre-operative 3D centerlines u. The single projected view consists of the 2D intra-operative centerlines v. The unknown 3D deformed shape δu is the displacement field induced by a potential brain shift. As in most registration methods, an initial set of n c correspondences c between the image points and the model points has to be established. We thus initialise our non-rigid registration with a rough manual rigid alignment. Once aligned, a vector c is built so that if a 2D point v i corresponds to a 3D point u j then c i = j for each point of the two sets. Assuming that the camera projection matrix P is known and remains constant, registering the 3D vessels to their corresponding 2D projections amounts to minimize the re-projection error Pu ci -v i for i ∈ n c . However, minimizing this error does not necessarily produces a correct 3D representation since many 3D shapes may have identical 2D projections. To overcome this issue, we add to the re-projection constraints the vessels' physical priors following the beams model introduced in Eq. 2. This will force the shape of the vascular tree to remain consistent in 3D while minimizing the reprojection error in 2D. This leads to the following force minimization expression:

argmin δu nc i f ci -κ Pu ci -v i 2 (3)
where κ is the stiffness coefficient that permits the image re-projection error to be translated to an image bending force. Note the subscript c i that denotes the correspondence indices between the two point sets. This minimization can be seen as enforcing 3D vessel centerlines to fit sightlines from the camera position to 2D vessels centerlines while maintaining a coherent 3D shape.

Mapping Tumors with Cortical Vessels Deformation

In order to propagate the surface deformation to tumors and other sub-cortical structures, we use a linear geometrical barycentric mapping function. We restrict the impact of the cortical vessel deformations to the immediate underlying structures (see Figure 6).

Formally speaking, if we denote the vector of vertices representing a 3D tumor by t, we can express each vertex t i using barycentric coordinates of facet vertices u, such that t i = 3 j=1 φ j (x i , y i , z i )u j , where φ(x, y, z) = a + bx + cy with (a, b, c) being the barycentric coordinates of the triangle composed of nodal points u j , with 1 ≤ j ≤ 3. This mapping is computed at rest and remains valid during the deformation.

Results

We tested our method on simulated data of a 3D synthetic human brain. The 3D vessel centerlines were extracted from the vessels mesh surface using a mean curvature skeletonization technique. The number of centerline nodes to build the wire-like model is set in a range between 12 and 30 nodes (only a subset of the vascular network has been considered). The stiffness matrix of the vessels is built with a Young's modulus set to 0.6 MPa and Poisson's ratio to 0.45 to simulate a quasi-elastic and incomprehensible behaviour where vessels thickness is set to 0.5 mm [START_REF] Ebrahimi | Mechanical properties of normal and diseased cerebrovascular system[END_REF]. We used the Sofa framework (www.sofa-framework.org) to build the 3D beam-elements model. We quantitatively evaluated our method by measuring the target registration error (TRE) on different locations. We considered four in-depth locations from 0 mm to 45 mm, 0 mm being the cortical surface. We simulated a brain shift to mimic a protrusion deformation that can occur due to brain swelling after the craniotomy. The deformations were of increasing amplitudes: Def1 = 1.7 mm, Def2 = 4 mm, Def3 = 6 mm and Def4 = 10 mm. We extracted 2D centerlines of each deformation by projecting the deformed 3D centerlines w.r.t a known virtual camera. We added Gaussian noise with a standard deviation of 5 mm and a 5 clustering decimation on the set of 2D points composing the centerlines. The results reported in Figure 3 suggest that using our method achieves small TRE ranging from 0.53 mm to 1.93 mm, on the cortical surface and the immediate subcortical structures (≤ 15 mm). Since the TRE depends on the amount of brain shift, we also quantify the brain-shift compensation by measuring the difference between the initial and the corrected error (normalized by the initial error). Brain shift compensation of up to 68.2 % is achieved, and at least a 24.6 % compensation is reported for the worst configuration. We can also observe that errors and compensation at the cortical surface are very close to the immediate sub-cortical location (≤ 15 mm). They increase when the targets are located deeper in the brain (≥ 30 mm) or the amount of deformation increases (≥ 6 mm). Fig. 4: Augmented Reality views after applying our method on 6 patient data sets retrospectively. The non-rigid registration provided a good estimation of deformation of the cortical surface, as viewed by overlays of the deformed 3D vessels on their 2D projections seen through the surgical microscope.

Our craniotomy dataset was composed of 1630 microscopic image patches with labels obtained through manual segmentation. The model is trained on mini-batches of size 32 for 200 epochs using Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with a learning rate of 0.001 that is decreased by a factor of 10 each 100 epochs. We used the Tensorflow framework (www.tensorflow.org) on an NVidia GeForce GTX 1070.

We tested our method on 6 patient data sets retrospectively. The cortical vessels, brain parenchyma, skull and tumors were segmented using 3D Slicer (www.sliced.org). The miscroscopic images were acquired with a Carl Zeiss surgical microscope. We only used the left image of the sterescopic camera.

The pre-operative vessels are first aligned manually on the image to obtain the rigid transform, then the 3D/2D non-rigid registration was performed, and used to drive deformation of the underlying tumors. Figure 4 shows the resulting Augmented Reality rendering and exhibits good mesh-to-image overlay. Figure 5 shows the measured deformation amplitude on each vessel after non-rigid registration. The deformations range from 0 to 3 mm and are non-uniformly distributed on the vascular trees which suggest a non-rigid deformation. Finally, we show in Figure 6 the ability of our method to rectify tumor positions through their mapping with the cortical vessels. We can clearly see that tumor positions are corrected according to the brain shift estimated from the cortical vessels registration.

Conclusion and Discussion

We proposed here a brain-shift aware Augmented Reality method for craniotomy guidance. Our method follows a 3D/2D non-rigid process between pre-operative and intra-operative cortical vessel centerlines using a single image instead of a stereo pair. Restricting our method to a single image makes it more acceptable in operating rooms but turns the registration process into an ill-posed problem.

To tackle this issue we proposed a force-based shape-from-template formulation that adds physical constraints to the classical reprojection minimization.

In addition, our pipeline takes advantage of recent advances in deep learning to automatically extract vessels from microscopic images. Our results show that low TRE can be obtained at cortical and sub-cortical levels and compensation for brain shift with up to 68% can be achieved. Our method is however sensitive to the outcome of the cortical vessels segmentation. A fragmented segmentation may lead to discontinuous centerlines and thus produce aberrant vascular trees that can affect the whole pipeline. One solution could be the use of the complete vessels segmentation in addition to centerlines to make advantage of the use of more featured such as radii, curves and bifurcations. Future work will consist of developing a learning-based method to perform the non-rigid registration without manual initialisation, with the aim of facilitating its usage by surgeons in clinical routines.

Fig. 1 :

 1 Fig.1: Problem formulation: we aim at recovering the deformed 3D vessels shape δu from its known reprojection in the image v, the known pre-operative 3D vessels at rest u and known rigid alignment c , satisfying physical and reprojective constraints P.

Fig. 2 :

 2 Fig. 2: Extracting 2D cortical brain vessels: A classical U-Net architecture is used to segment the microscopic images. The input images (a) or (d) are segmented to obtain a binary image (b) or (e). The images (c) and (f) consist of overlays.

Fig. 3 :

 3 Fig. 3: Quantitative evaluation on simulated data: Charts (a) and (b) show the percentage of compensation and target registration error respectively w.r.t to depth and brain shift amplitude. The depth axis correspond to the position of the target in the brain where 0mm represent the cortical surface and 45mm the deep brain.

Fig. 5 :

 5 Fig. 5: Color-coded estimated deformation on the vessels.

Fig. 6 :

 6 Fig. 6: Mapping tumors with cortical vessels: In (a) and (c), a visualization of the mapping mechanism. Blue lines represent the attachments. In (b) and (d), an Augmented Reality visualization of the tumors with the estimate displacement induced from the mapping w.r.t to initial positions (dotted blue lines).

Acknowledgements

The authors were supported by the following funding bodies and grants: NIH: R01 EB027134-01, NIH: R01 NS049251 and BWH Radiology Department Research Pilot Grant Award.