Viet Hoang Man

Xibing He

Phuong H Nguyen

Celeste Sagui

Christopher Roland

Xiang-Qun Xie

Junmei Wang
email: junmei.wang@pitt.edu

Unpolarized laser method for infrared spectrum calculation of amide I C O bonds in

Keywords: IR spectra, Amide I, radiation absorption, external electric field, MD simulation

  L'archive ouverte pluridisciplinaire

INTRODUCTION

Infrared (IR) spectroscopy, which is based on the absorption of IR radiation, is one of the most popular methods in molecular structure determination. [START_REF] Barth | Infrared spectroscopy of proteins[END_REF] It has been widely used to investigate protein structures in protein folding, unfolding and misfolding. [START_REF] Kauffmann | Lifetimes of intermediates in the β-sheet to α-helix transitions of β-lactoglobulin by using a diffusional IR mixer[END_REF][START_REF] Fabian | Two-dimensional IR correlation spectroscopy: Sequential events in the unfolding process of the lambda Cro-V55C repressor protein[END_REF][START_REF] Reinstädler | New structural insights into the refolding of ribonuclease T1 as seen by time-resolved Fourier-transform infrared spectroscopy[END_REF][START_REF] Schultz | Illuminating folding intermediates[END_REF][START_REF] Troullier | Transient non-native secondary structures during the refolding of -lactalbumin detected by infrared spectroscopy[END_REF][START_REF] Ramajo | The α-helix folds more rapidly at the Cterminus than at the N-terminus[END_REF] Among the IR normal modes of the protein IR spectra, the amide I mode, which arises mainly from the C=O stretching vibration, is most commonly used for protein secondary structure analysis due to its strong signal. Additionally, the amide I vibration is hardly affected by the nature of the side-chain, but it is sensitive to hydrogen bonding, dipole-dipole interactions, and the conformation of the backbone. [START_REF] Krimm | Vibrational Spectroscopy and Conformation of Peptides, Polypeptides, and Proteins[END_REF][START_REF] Henry | Infrared Spectroscopy of Biomolecules[END_REF] The amide I IR spectra of a protein for a given protein structure can reveal the percentages of the basic secondary structural contents including α-helix, β-sheet, turns and disordered loops in the protein structure. [START_REF] Byler | Examination of the secondary structure of proteins by deconvolved FTIR spectra[END_REF][START_REF] Barth | What vibrations tell us about proteins[END_REF] However, it cannot directly show the protein tertiary structure which may be modeled from the data of X-ray or NMR experiments, or by pure theoretical modeling. In this context, one can first calculate the theoretical IR spectra of a protein tertiary structure model, and then investigate how well the calculated IR spectra are consistent with the experimental ones. In this way, the quality of the protein model can be assessed.

Moreover, the combination of theoretical and experimental amide I IR spectra can reveal the transition-state structures in a protein dynamics process. [START_REF] Reppert | Computational Amide I 2D IR Spectroscopy as a Probe of Protein Structure and Dynamics[END_REF][START_REF] Davis | A quantitative connection of experimental and simulated folding landscapes by vibrational spectroscopy[END_REF] In the last few decades, several theoretical-computational approaches have been proposed to reproduce the amide I IR spectra of well-defined molecules. [START_REF] Torii | Model calculations on the amide-I infrared bands of globular proteins[END_REF][START_REF] Lee | Ab initio-based vibrational analysis of α-poly(L-alanine)[END_REF][START_REF] Brauner | A quantitative reconstruction of the amide I contour in the IR spectra of globular proteins: from structure to spectrum[END_REF][START_REF] La Cour Jansen | Modeling the amide I bands of small peptides[END_REF][START_REF] Smith | Amide I two-dimensional infrared spectroscopy of β-hairpin peptides[END_REF][START_REF] Ham | Amide I modes of tripeptides: Hessian matrix reconstruction and isotope effects[END_REF][START_REF] Hahn | Simulation studies of amide I IR absorption and twodimensional IR spectra of β hairpins in liquid water[END_REF][START_REF] Choi | Computational spectroscopy of ubiquitin: comparison between theory and experiments[END_REF][START_REF] Wang | Development and Validation of Transferable Amide I Vibrational Frequency Maps for Peptides[END_REF][START_REF] Polzi | Analysis of Infrared Spectra of β-Hairpin Peptides As Derived from Molecular Dynamics Simulations[END_REF][START_REF] Xue | Perturbation approach for computing infrared spectra of the local mode of probe molecules[END_REF][START_REF] Olson | Enhanced vibrational solvatochromism and spectral diffusion by electron rich substituents on small molecule silanes[END_REF][START_REF] Yin | Active-Site Heterogeneity of Lactate Dehydrogenase[END_REF] These approaches are typically based on mixed quantum-classical strategies and can be divided into two categories. In the first category, the vibration frequencies and eigenstates for single amides are determined by employing quantum mechanical calculations and then calibrated for the whole protein. [START_REF] Torii | Model calculations on the amide-I infrared bands of globular proteins[END_REF][START_REF] Lee | Ab initio-based vibrational analysis of α-poly(L-alanine)[END_REF][START_REF] Brauner | A quantitative reconstruction of the amide I contour in the IR spectra of globular proteins: from structure to spectrum[END_REF][START_REF] La Cour Jansen | Modeling the amide I bands of small peptides[END_REF][START_REF] Smith | Amide I two-dimensional infrared spectroscopy of β-hairpin peptides[END_REF] In the second category the Hessian matrix for the whole isolated protein is calculated, and then the local, single-residue-based, vibrational frequencies are obtained via a matrix reconstruction. [START_REF] Ham | Amide I modes of tripeptides: Hessian matrix reconstruction and isotope effects[END_REF][START_REF] Hahn | Simulation studies of amide I IR absorption and twodimensional IR spectra of β hairpins in liquid water[END_REF][START_REF] Choi | Computational spectroscopy of ubiquitin: comparison between theory and experiments[END_REF] In these approaches, diagonalization of the whole Hamiltonian matrix or localized amide Hamiltonian matrices provides the vibrational eigenstates and energy of the whole protein for a given protein conformation. These approaches have been successfully applied to investigate different proteins, including β-hairpin peptides, ubiquitin, and other globular proteins, as well as to characterize the IR spectral features of individual secondary-structure elements. [START_REF] Torii | Model calculations on the amide-I infrared bands of globular proteins[END_REF][START_REF] Brauner | A quantitative reconstruction of the amide I contour in the IR spectra of globular proteins: from structure to spectrum[END_REF][START_REF] Smith | Amide I two-dimensional infrared spectroscopy of β-hairpin peptides[END_REF][START_REF] Ganim | Amide I two-dimensional infrared spectroscopy of proteins[END_REF] Conversely, given an IR spectrum, we can deduce the secondary-structure elements that contributed to it.

In this study, we simulated an "unpolarized" laser via external electric field (EEF) to probe the infrared spectrum of C=O bonds in different protein structures by using all-atom nonequilibrium molecular dynamics (MD) simulations. In our method, the IR radiation is mimicked by a time dependent EEF which has a shape of laser pulse and IR frequencies. The EEF is applied along the directions of the individual C=O bonds to get maximize EEF's effect.

The basic principle of the method is that the EEF causes C=O stretching vibration and subsequentially leads to a system energy change which represents the radiation absorption. This application of an EEF along covalent bonds allows us to simulate single-wavelength radiation.

Note that polarized lasers are popular in experiments as well as applications for IR spectroscopy.

In an experimental IR spectrum measurement, each protein in the solution will only experience an electric field in a given direction, but the signal is averaged over the ensemble. Therefore, to model an experiment with a polarized laser, we need to use many same proteins with different orientations in a simulation system and/or multiple scanning with different angles of the EEF. This will result in considerable computational cost. On the other hand, the simulated "unpolarized" laser, where all excitation orientations are sampled, allows us to model an experimental IR measurement with a single protein in a simulation system and without the consideration of the EEF angles. This method has been implemented in the AMBER software package [START_REF] Case | AMBER 19[END_REF] . Using the unpolarized radiation, we computed the amide I IR spectrum of the C=O bonds in representative protein structures including β-sheet and α-helix. The optimization of the important input parameters for the IR spectrum calculation has been investigated. Finally, we discuss the application of the laser radiation for force field development and molecular mechanism studies for radiation absorption.

METHODS AND MATERIALS

Protein structural models and equilibrium MD simulation

The protein structures selected to calculate C=O bond infrared spectra include amyloid antiparallel β-sheet (PDB ID: 2KIB), [START_REF] Nielsen | Unique Identification of Supramolecular Structures in Amyloid Fibrils by Solid-State NMR Spectroscopy[END_REF] amyloid parallel β-sheet (PDB ID: 1YJP), [START_REF] Nelson | Structure of the cross-β spine of amyloid-like fibrils[END_REF] β-hairpin (Tryptophan Zipper 2, PDB ID: 1LE1), [START_REF] Cochran | Tryptophan zippers: Stable, monomeric β-hairpins[END_REF] α-helix (PDB ID: 6N8C), [START_REF] Kotler | Probing initial transient oligomerization events facilitating Huntingtin fibril nucleation at atomic resolution by relaxation-based NMR[END_REF] and the structure of the transmembrane region of the M2 protein H + channel (PDB ID: 1MP6). [START_REF] Wang | Structure of the transmembrane region of the M2 protein H+ channel[END_REF] The representative protein secondary structures, anti-parallel β-sheet (PDB ID: 2KIB), parallel β-sheet (PDB ID: 1YJP) and α-helix (PDB ID: 6N8C), are shown in Figure 1b. For each structure, the protein was placed at the center of an octahedron box with a minimum distance of 10 Å between box borders and any atoms of protein. The protein was solvated with explicit water using the TIP3P water model. [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF] The ff14SB force field [START_REF] Maier | ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB[END_REF] was used to describe proteins. The Na + /Cl -ions were properly added to neutralize the systems. Detailed information of the simulation systems is provided in the Supporting Information (Table S1). Each MD system underwent the following steps. First, the steepest descent minimization followed by a conjugate gradient minimization with the peptide atoms restrained at their initial positions, then the steepest descent minimization without restraint followed by a conjugate gradient minimization were carried out. The minimization step was followed by a short MD simulation under constant volume while the system was heated from 0 K to 300 K with weak restraints on the protein atoms. Next, Langevin dynamics at constant temperature (300 K) and constant pressure (1 atm) were run for 100 ps, after which the density of the system was found to be stable around 1.0 g/cm [START_REF] Fabian | Two-dimensional IR correlation spectroscopy: Sequential events in the unfolding process of the lambda Cro-V55C repressor protein[END_REF] . Finally, in the sampling phase, NPT simulations at 300 K were performed using the leap-frog algorithm with a time step of 2 fs and a duration time of 10 ns. The long-range Coulomb interaction was accounted by means of the Particle-Mesh Ewald (PME) method [START_REF] Essmann | A smooth particle mesh Ewald method[END_REF] with a 10 Å cutoff. The van der Waals interactions were calculated by means of a 10 Å atom-based nonbonded lists, with continuous corrections applied to the long-range parts. The constant pressure simulations were carried out at 1 atm via the Berendsen barostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] with the pressure relaxation time τ p of 3.0 ps. The temperature was regulated using Langevin dynamics with a collision frequency of 1 ps -1 . The SHAKE algorithm [START_REF] Forester | rattle, and roll: Efficient constraint algorithms for linked rigid bodies[END_REF] was applied to all bonds involving hydrogen atoms. The periodic boundary condition was applied to all MD simulations which were performed using the pmemd.cuda module implemented in the AMBER 19 package. [START_REF] Case | AMBER 19[END_REF]

Simulated laser applied on the amide bond (C=O)

An "unpolarized" laser pulse was simulated by a time dependent electric field as shown in Equation 1.

(1)

Here, (in V/Å) represents the amplitude of the electric field, (ps) is the pulse temporal width, t (ps) is the time, (ps) is the time when the pulse is the maximum (E 0 ), c (cm/ps) is the speed of light, and (cm -1 ), the space frequency is the inverse of the wavelength, as such 2π is the wavenumber. This technique has been implemented in the GROMACS simulation package, [START_REF] Lindahl | GROMACS 3.0: a package for molecular simulation and trajectory analysis[END_REF] but not in the AMBER package [START_REF] Case | AMBER 19[END_REF] yet. We have used this mimic laser for laser-induced simulations to study the dissociation of amyloid fibrils, poliovirus, peptide nanotubes and cellulose fibers and fast DNA melting. [START_REF] Man | Picosecond Dissociation of Amyloid Fibrils with Infrared Laser: A Nonequilibrium Simulation Study[END_REF][START_REF] Man | Picosecond infrared laser-induced all-atom nonequilibrium molecular dynamics simulation of dissociation of viruses[END_REF][START_REF] Man | Picosecond Melting of Peptide Nanotubes with Infrared Laser: A Nonequilibrium Simulation Study[END_REF][START_REF] Man | Comparative melting and healing of B-DNA and Z-DNA by an infrared laser pulse[END_REF][START_REF] Domin | Breaking Down Cellulose Fibrils With a Mid-infrared Laser[END_REF] The input parameters for a laser-induced simulation include , , and , in which, and control the laser power. In these laser simulations, the mimic laser is applied on a molecule such that the harmonic electric field excites the vibrational transitions of ions or partially charged atoms with the corresponding laser frequency, . In the case of the C=O amide bond, C and O atoms have positive and negative partial charges, respectively. Therefore, their vibration transitions are opposite to each other, resulting in bond stretch vibration. If the electric field frequency, , equals the natural frequency of the bond, the bond vibration will be resonant, and the energy of the simulated system will increase. In other words, one can state that the amide bond absorbs the laser energy at the specific bond stretch frequency, namely amide I frequency. In experiment, although the polarized laser is used, there are many molecules in the sample, and the measurement is repeated for each sample with different angles. Then, the spectrum is basically averaged over all orientations of molecules. In this context, we can think that the laser used in experiment is "unpolarized". On the other hand, in MD simulations, the sample usually consists of a single protein. Therefore, if the mimic laser behaves like a polarized light beam for which the electric field vectors only vibrate in a specific plane, it shows maximum effect on the vibration of a bond only when the electric field vector is applied along with the bond direction, and it has no or very little effect on the bond vibration when its electric field vector is perpendicular to the bond direction. The mimic laser can be used to measure the infrared spectra of amide I in some special peptide/protein structures such as βsheet and single helices in which all C=O bonds are mostly in the same direction. In those cases, we can choose the direction of the electric field vector as parallel to the C=O bond to maximize the laser effect. However, in most protein structures, all C=O bonds are not ordered and take random directions. Thus, it is impossible to mimic the effect of applying a polarized laser on an ensemble of protein molecules in experiment by using a single protein in macromolecule in periodic MD simulations.

To overcome the above limitation, we applied a previously unused strategy to mimic the effect of applying an unpolarized laser to a protein: for a given C=O bond, the magnitude of the electric field is described by Equation 1, but the direction of the electric field vector is adapted and always aligned to the direction of the C=O bond. With this new strategy, the forces acting on the C and O atoms due to the electric field, and , are now calculated by the Equations 2-3:

(

Where, is electric field magnitude described by Equation 1. and are partial charges of C and O atoms, respectively. The is the unit vector which points from C atom to O atom. The parameters for controlling the laser now include the C=O bond information, , , and . By scanning some parameters especially , we are able to study the infrared absorption of a covalent bond mimicking the effect of applying a polarized laser beam on a protein sample in experiment.

Laser-induced simulations for the infrared spectra investigation

To investigate the IR spectra of amide bonds in the protein structures, we have implemented aforementioned technique into the AMBER software package. For each considered protein structure, we performed the laser-induced simulations for 201 parameter sets with the laser frequencies, , varied from 1550 cm -1 to 1750 cm -1 at a step of 1 cm -1 and the other laser parameters being fixed at values of = 0.18 V/Å, = 3 ps and = 1 ps. The total simulation time is 10 ps. For each parameter set/frequency, ten trajectories were carried out with ten different starting conformations taken from the 10-ns equilibrium simulation. We also investigated ten different simulated time steps varying from 0.01 fs to 2 fs for each trajectory.

Thus, there are about 20,100 laser-induced simulations which were carried out for each protein model. Additionally, we carried out the laser-induced simulation to estimate the effect of C=O bond stretching force constant on the IR spectra (for the simulation time step of 0.1 fs only).

Overall, the total number of laser-induced simulations performed in this study is more than 80,000.

Data analysis

The time-dependent fluctuation of a C=O bond was calculated by the following equation:

, where is time evolution of the i th C=O bond and 1.23 Å is its equilibrium value. The time-dependent fluctuation of the C=O bond in a protein is the ensemble average, , which was averaged over all trajectories and all C=O bonds in the protein.

To estimate the IR radiation absorption of a protein system, we calculate the time-dependent difference between the laser-induced and the equilibrium energies: S . Here, is the energy of the system in a laser-induced simulation at a given time, and is the averaged equilibrium energy at the same temperature as in the laser-induced simulations.

RESULTS AND DISCUSSIONS The IR radiation absorption of C=O bonds

The basic principle of IR spectroscopy is that when IR radiation passes through a sample, the absorption of IR radiation by the sample produces a change of the radiation intensity, and the change is then detected and translated into the absorption intensity of the sample. In our laserinduced simulation, the laser mimicked by the time-dependent electric field (Figure 1a) induces the vibration of the targeted bonds, resulting in a change of the system energy. Therefore, the IR absorption can be described by the difference between the energies of the system before and after the laser pulse applied. Figure 2 shows the time dependence of the energy differences and the fluctuation of the C=O bonds in the different protein systems during the laser-induced simulations for the given laser frequencies. The data shows that the energy difference and the bond fluctuation strongly depend on the laser frequency, and they are highly correlated. For example, for laser frequencies of 1630 cm -1 and 1720 cm -1 , the energy of all considered systems did not change (or changed very little). Their C=O fluctuation is very small, suggesting that the C=O bonds did not absorb the laser radiation energy at those frequencies. In contrast, for laser frequencies of 1676 cm -1 , 1679 cm -1 and 1691 cm -1 , the C=O fluctuation is clearly shown. The energy difference was significantly changed during the simulation time. The time dependence of the energy difference can be divided into three states. In the first state, from the beginning of the simulation to 2 ps, the energy mostly did not change. The second state, from 2 ps to 4 ps, the energy difference dramatically increased. In the last state, the energy difference slightly increased from 4 ps and reached a plateau value at about 5 ps. Therefore, for each laser frequency, the energy difference at the last 10 th ps of the laser-induced simulation will be used as the protein absorption energy, which is equal to the absorption intensity in IR spectra.

Important input parameters for the laser-induced simulation

Since the mimic laser absorption in the simulation was caused by the C=O vibration, the temperature coupling parameters, which directly affect the instantaneous velocities of the atoms, are important to a laser-induced simulation. In principle, when a system is coupled to the thermostat, the velocities of all atoms in the systems will be rescaled to get the system temperature close with the desired value. This implies that the effect of laser on the atom vibrations will diminish when the temperature coupling is applied. Therefore, in an ideal laserinduced simulation, the system is only coupled to the thermostat at the beginning of the simulation for its initial thermostat temperature, and the temperature coupling should not be applied during the rest of the simulation time. To satisfy the temperature coupling condition, in our laser-induced simulations, the flag ntt, which is used for the temperature regulation, was set to 3 and with gamma_ln being set to 0.00001. This setting will rescale velocities of the atoms every 100,000 ps, which is much longer than the 10 ps of a laser-induced simulation.

In the laser-induced simulations, the external forces due to the electric field applied on the atoms of the C=O bonds are calculated by Equations 2 and 3. Because the partial charge of an atom is constant, the time dependence of the external electric force is similar to that of the external electric field described by Equation 1. In any MD simulation, the force acting on an atom and the velocity of the atom is calculated and updated for every simulated time step, namely dt. Therefore, the external electric force contributed by laser is discrete. If the laser period, , is too small, the ratio T/dt (resolution) will be significant small, the effect of the discreteness will become obvious, and the external electric field will be very different from the one truly applied in the simulations. The effect of dt on a mimic laser pulse, which has the frequency of 6000 cm -1 (corresponding to the period T of 5.56 fs) was shown in Figure 3.

It shows that when dt = 0.02 fs, the resolution T/dt is large enough to reproduce the periodicity of the laser, while the resolution will be too low for the laser expression if dt = 2 fs. In general, the higher the resolution is, the better the time-dependent pattern of the laser is reproduced, and the unpolarized laser effect better achieved. Although applying a small dt can increase the resolution, it also considerably slows the simulation. Thus, it is important to determine reasonable values of dt for laser-induced simulations at a given range of laser frequencies. In this study, the scanning frequencies were in the range of 1550 cm -1 to 1750 cm -1 . We have measured to low frequencies, then the peak is unchanged as dt is smaller than 0.4 fs (Figure 4d). In other words, a dt value smaller than 0.4 fs is good enough for studying the amide I spectra via laserinduced simulations. Later, we will discuss the results of simulations using dt =0.1 fs.

Other tunable parameters include and , which control the energy of a laser pulse at a given frequency. If the laser energy is too high, the structure of the protein may dramatically change, discrepancies may also occur because the calculated spectrum is not for the true protein structure at all. If the laser energy is too low, the absorbed peak will be very low, which may not be distinguishable from the thermal noise. As the demonstration in Figure 5, when = 0.5

V/nm and = 1 ps, the absorption energy is very small, while if = 3.0 V/nm and = 3 ps the protein structure significantly changed in comparison with the initial conformation. In this study, we used the laser parameters, = 1.8 V/nm and = 1 ps. It provides a good observation of the IR radiation absorption but has neglectable effect on the structure of protein models.

Amide I spectra of proteins

Amide I spectra of secondary structural elements

One of the most important applications of a protein amide I IR spectra is to determine protein secondary structures. In this application, the correct identification of the absorbed peak positions, which correspond to the protein secondary structure elements, is crucial. Goormaghtigh et al.

have collected and evaluated the assignment of amide I band positions to protein secondary structure based on experimental data and assignments in various publications. [START_REF] Goormaghtigh | Determination of Soluble and Membrane Protein Structure by Fourier Transform Infrared Spectroscopy III. Secondary Structures[END_REF] Their study suggested that two IR regions, 1623-1641 cm -1 and 1674-1695 cm -1 , are assigned for β-sheet structures, and a region, 1648-1657 cm -1 , is assigned for α-helix structures. The average positions of those assignments are 1633 cm -1 and 1684 cm -1 for β-sheets, and 1654 cm -1 for α-helixes.

Those assignments have been widely applied in analyzing the Fourier-transform infrared (FTIR) spectroscopy. [START_REF] Lindahl | GROMACS 3.0: a package for molecular simulation and trajectory analysis[END_REF] In this study, we applied our specific method to calculate amide I spectra of some protein models which contain only one type of secondary structure elements, and we compared our calculated spectra with the average experimental data.

Figure 6 shows the amide I IR spectra of the antiparallel β-sheet, parallel β-sheet and α-helix protein structures which were calculated by nonequilibrium MD simulations with the external electric field mimicking a laser pulse, with = 1.8 V/nm, = 1 ps and dt = 0.1 fs. For the βsheet structures, our calculated spectra contain only one peak. The absorbed peak of the parallel β-sheet spectra is located at 1675 cm -1 , and the peak for the antiparallel case slightly shifted to a low frequency, and 5-9 cm -1 lower than the experimental high frequency, respectively. [START_REF] Goormaghtigh | Determination of Soluble and Membrane Protein Structure by Fourier Transform Infrared Spectroscopy III. Secondary Structures[END_REF] Many previous studies agreed that the spectra of both the β-sheet structures contains a main peak with strong intensity at low frequency region (1617-1641 cm -1 ), and the peak with much lower intensity at the high frequency region (above 1680 cm -1 ) is a feature of antiparallel β-sheet structure only. [START_REF] Zandomeneghi | FTIR reveals structural differences between native -sheet proteins and amyloid fibrils[END_REF][START_REF] Fu | In Situ Misfolding of Human Islet Amyloid Polypeptide at Interfaces Probed by Vibrational Sum Frequency Generation[END_REF][START_REF] Kubelka | Differentiation of β-Sheet-Forming Structures: Ab Initio-Based Simulations of IR Absorption and Vibrational CD for Model Peptide and Protein βSheets[END_REF][START_REF] Zou | Parallel β-Sheet Fibril and Antiparallel β-Sheet Oligomer: New Insights into Amyloid Formation of Hen Egg White Lysozyme under Heat and Acidic Condition from FTIR Spectroscopy[END_REF] Thus, our observed peaks should belong to the low frequency region, but they were shifted 40 cm -1 toward the higher frequency direction in comparison with the average experimental value.

We did not observe the low intensity peak at high frequency region in our calculated spectra. For the α-helix structure, we observed a peak at 1692 cm -1 , which is again higher than the experimental observation about 40 cm -1 . Despite of the significant shift of our calculated spectrum peaks compared to the experimental values, the difference between α-helix and β-sheet peaks in our calculation (about 17 cm -1 ) is comparable with the experimental finding (about 20 cm -1 ). Additionally, our result, which showed a slightly difference between the spectra of parallel and antiparallel β-sheet structures, is consistent with a previous theoretical calculation. [START_REF] Roeters | Evidence for Intramolecular Antiparallel Beta-Sheet Structure in Alpha-Synuclein Fibrils from a Combination of Two-Dimensional Infrared Spectroscopy and Atomic Force Microscopy[END_REF]

Correction on the C=O bond force-constant of the FF14SB force field

Although our IR spectra were calculated by targeting only the vibration of the C=O bond, we still observed the distinguished peaks for different protein structures as experiments and other theoretical calculations did. However, the peak positions in our spectrum calculations were different from the ones assigned by experimental data. We have discussed how the IR spectra are affected by the laser and simulation parameters. We showed that the peak positions only depend on the simulation time step dt when dt is larger than 0.4 fs, and this dependency disappear when dt is smaller than 0.4 fs. Thus, we wondered whether the difference was caused by force field parameters. In order to investigate the factors that affect the absorbed frequencies of the C=O bonds in different protein conformations, let us consider a simple two-atom oscillator, which has the natural frequency ( ) given by the equation:

(4) (5)
where, the natural frequency is in cm -1 , c is speed of light, is the force constant between the two atoms, is the reduced mass, is mass of the first atom, and is mass of the second atom. Therefore, two factors influencing the frequency are the force constant and the mass of the vibrating atoms. The frequency rises when the force constant increases and/or the reduced mass decreases, and vice versa. Applying this oscillator to the case of our spectrum calculation with all-atom MD simulation, is the force constant of C=O bond, which was defined by the force field and did not change during the simulations, and m 1 and m 2 are effective masses of C and O atoms. Those effective masses mostly depend on the masses of the atoms, intermolecular interactions, intramolecular interaction and conformation of the protein, resulting in the dependence of IR absorbed peak on the protein structure. Note that, like the force constant of a bond, the mass of an atom was also constant and defined by force field. In the AMBER FF14SB [START_REF] Maier | ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB[END_REF] we used, the atomic mass was experimentally measured, while the force constant was usually calculated using quantum chemical methods. Therefore, we only consider the impact of the bond force constant on the protein spectra. As seen in Figure 7, the absorbed peaks of both antiparallel β-sheet and α-helix spectra shifted to lower frequency when the C=O bond force constant was reduced. Interestingly, when =515 kcal/mol/Å 2 the absorbed peaks of the βsheet and α-helix spectra are 1633 cm -1 and 1654 cm -1 , which are the same as the experimental averaged values assigned to the β-sheet and α-helix secondary structure. In other words, with the C=O bond stretching force constant, being set to 515 kcal/mol/Å 2 , not the original one (570 kcal/mol/Å 2 ) in the AMBER ff14SB force field, our method can exactly reproduce the amide I IR spectra of β-sheet and α-helix protein structures. Therefore, our method can also be applied to adjust the bond force constant parameters of a force field, which may not be ideal for all-atom MD simulations of biomolecules as those force field parameters were usually developed through the quantum mechanical calculation of prototype molecules in vacuum. Note that most of AMBER biomolecular force fields share the same bond force-constants of the bonds in proteins.

However, other force fields such as GROMOS [START_REF] Hermans | A consistent empirical potential for water-protein interactions[END_REF] or OPLS-AA [START_REF] Kaminski | Evaluation and Reparameterization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides[END_REF] do not share the bond-force constants with AMBER force fields. Therefore, the correction of the bond-force constants based on IR spectra should be done for an individual force field or for a force field family.

Application of unpolarized laser all-atom MD simulations to calculate amide I spectra for representative protein secondary structures.

In this section, we applied the mimic laser to calculate the IR spectra of additional proteins and compared the calculated spectra to the experimental data. It is noted that the force constant of C=O bond stretching parameter, , was set to 515 kcal/mol/Å 2 . The first system considered was Tryptophan Zipper2 (Trpzip2), a β-hairpin structure. The experimental amide I spectrum of Trpzip2 was produced by Smith and Tokmakoff. [START_REF] Smith | Amide I two-dimensional infrared spectroscopy of β-hairpin peptides[END_REF] Our calculated spectrum agreed with the experiment in terms of the peak location that was around 1636 cm -1 (Figure 8). However, the calculated spectrum peak width was narrower than the experimental one, probably due to insufficient conformational sampling.

The second protein system we considered was the three peptides encompassing the transmembrane domain of Influenza A virus M2 protein H + channel (Ser22-Leu46, SSDPLVVAASIIGILHLILWILDRL). The first peptide was not labeled (wild type peptide), the second one contains a 13 C- [START_REF] Brauner | A quantitative reconstruction of the amide I contour in the IR spectra of globular proteins: from structure to spectrum[END_REF] O isotopic labeled alanine residue at position 29 ( 13 C- [START_REF] Brauner | A quantitative reconstruction of the amide I contour in the IR spectra of globular proteins: from structure to spectrum[END_REF] O A29 peptide), and the third peptide contains a 13 C- [START_REF] Smith | Amide I two-dimensional infrared spectroscopy of β-hairpin peptides[END_REF] O isotopic labeled glycine residue at position 34 ( 13 C-18 O G34 peptide). The amide spectra of those peptides were experimentally examined by Torres et al. for studying the vital role of the isotopic labels in relaying site-specific secondary structure and orientational information of proteins. [START_REF] Torres | Site-specific examination of secondary structure and orientation determination in membrane proteins: The peptidic 13 C-18 O group as a novel infrared probe[END_REF] To calculate the peptide spectra, the initial structure of the wild-type peptide was taken from PDB databank (PDB ID: 1MP6), while the isotopic labeled residues were obtained by properly changing the mass of the related atom in the AMBER topological file of the wild-type peptide. All systems of the peptides were created and underwent the simulated process as described in the method section. As illustrated in Figure 9, the calculated spectra were close to experimental ones. The calculated spectrum peaks, which related to [START_REF] Reppert | Computational Amide I 2D IR Spectroscopy as a Probe of Protein Structure and Dynamics[END_REF] C- [START_REF] Brauner | A quantitative reconstruction of the amide I contour in the IR spectra of globular proteins: from structure to spectrum[END_REF] O bond absorptance, was at 1654 cm -1 , differing from experiment (1659 cm -1 ) by only 5 cm -1 . The experimental weak peak at 1619 cm -1 associated with the isotopic [START_REF] Davis | A quantitative connection of experimental and simulated folding landscapes by vibrational spectroscopy[END_REF] C- [START_REF] Brauner | A quantitative reconstruction of the amide I contour in the IR spectra of globular proteins: from structure to spectrum[END_REF] O labeled A29 was not clearly observed in our calculation, while the experimental weak peak at 1595 cm -1 in 13 C-18 O G34 peptide system was reproduced in our theoretical spectrum. This redshifted peak can be qualitatively explained by the increasing of masses of labeled atoms (Equation .2). In other words, our method was able to describe the experimental red-shifted spectra caused by isotopic labeling if the effect is significant, such as the [START_REF] Davis | A quantitative connection of experimental and simulated folding landscapes by vibrational spectroscopy[END_REF] C-18 O isotopic label.

Last, we evaluated our method using a much larger protein system, the receptor binding domain of the spike protein of SARS-CoV2. We calculated the spectrum maps of all the residues of the protein in two different states, the apo state and the binding state in complex with the ACE2 protein as shown in Figure 10. The initial structure of the RBD-SARS-CoV2/ACE2 complex was taken from PDB databank (PDB ID: 6M0J). [START_REF] Lan | Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor[END_REF] Two MD systems, the RBD-SARS-CoV2/ACE2 complex and the apo-structure of RBD-SARS-CoV2 (extracted from the complex), were prepared and simulated following the same protocols as described in the Method section.

The residue-spectrum maps of RBD-SARS2 protein in the binding and unbound state were constructed and compared (Figure 10). The difference of the maps reveals the residues which are possibly keys to the RBD-SARS2 and ACE2 binding. Interestingly, many residues, such as 73D, 74E, 114G, 143A, 144G, 161Q, 162S and 165F, which have the largest differences in vibrational frequencies between the two states are located in the Spike and ACE2 binding interface. Those residues could be hotspot residues for the protein-protein interactions.

Interestingly, many of those residues are also recognized as hotspot residues through free energy decomposition analysis including 143A and 161Q. Still the neighbors of some residues like 144G, 162S, 165F are hotspot residues in free energy decomposition analysis. 56

Limits of unpolarized laser all-atom MD simulations in calculating IR spectra

A number of methods for IR spectra calculation, remarkably vibrational frequency maps, have been recently developed, and most of them are based on the mixed quantum-classical treatment. [START_REF] La Cour Jansen | Modeling the amide I bands of small peptides[END_REF][START_REF] Smith | Amide I two-dimensional infrared spectroscopy of β-hairpin peptides[END_REF][START_REF] Ham | Amide I modes of tripeptides: Hessian matrix reconstruction and isotope effects[END_REF][START_REF] Hahn | Simulation studies of amide I IR absorption and twodimensional IR spectra of β hairpins in liquid water[END_REF][START_REF] Choi | Computational spectroscopy of ubiquitin: comparison between theory and experiments[END_REF][START_REF] Wang | Development and Validation of Transferable Amide I Vibrational Frequency Maps for Peptides[END_REF][START_REF] Polzi | Analysis of Infrared Spectra of β-Hairpin Peptides As Derived from Molecular Dynamics Simulations[END_REF][START_REF] Xue | Perturbation approach for computing infrared spectra of the local mode of probe molecules[END_REF][START_REF] Olson | Enhanced vibrational solvatochromism and spectral diffusion by electron rich substituents on small molecule silanes[END_REF] Those state-of-the-art methods typically utilize biomolecular simulations to calculate IR spectra and the theoretical IR spectra can quantitatively reproduce the experimental ones. [START_REF] Reppert | Computational Amide I 2D IR Spectroscopy as a Probe of Protein Structure and Dynamics[END_REF][START_REF] Wang | Development and Validation of Transferable Amide I Vibrational Frequency Maps for Peptides[END_REF][START_REF] Olson | Enhanced vibrational solvatochromism and spectral diffusion by electron rich substituents on small molecule silanes[END_REF][START_REF] Yin | Active-Site Heterogeneity of Lactate Dehydrogenase[END_REF] However, the mixed quantum-classical approaches are usually complicated and their broad applications may be limited. To our best knowledge, the vibrational frequency maps take into account only the second nearest coupling between amide groups, and the contribution of the solvent is usually calculated based on classical MD simulation trajectories of small dipeptides or tripeptides, using the electric field or potential methods (e. g. Skinner's map). [START_REF] Wang | Development and Validation of Transferable Amide I Vibrational Frequency Maps for Peptides[END_REF] Therefore, these vibrational frequency maps produce very good IR spectra of small peptides, but for large biomolecules such as proteins and amyloid fibrils, the accuracy of these maps is questionable. This situation indicates that the simulation of IR spectra of biomolecules is highly challenging. In addition, the vibrational frequency maps have been parameterized only for a few biomolecular building blocks, including amino acids, [START_REF] La Cour Jansen | Modeling the amide I bands of small peptides[END_REF][START_REF] Smith | Amide I two-dimensional infrared spectroscopy of β-hairpin peptides[END_REF][START_REF] Ham | Amide I modes of tripeptides: Hessian matrix reconstruction and isotope effects[END_REF][START_REF] Hahn | Simulation studies of amide I IR absorption and twodimensional IR spectra of β hairpins in liquid water[END_REF][START_REF] Choi | Computational spectroscopy of ubiquitin: comparison between theory and experiments[END_REF][START_REF] Wang | Development and Validation of Transferable Amide I Vibrational Frequency Maps for Peptides[END_REF][START_REF] Polzi | Analysis of Infrared Spectra of β-Hairpin Peptides As Derived from Molecular Dynamics Simulations[END_REF][START_REF] Xue | Perturbation approach for computing infrared spectra of the local mode of probe molecules[END_REF][START_REF] Olson | Enhanced vibrational solvatochromism and spectral diffusion by electron rich substituents on small molecule silanes[END_REF] and nucleobases. [START_REF] Jiang | Development of Vibrational Frequency Maps for Nucleobases[END_REF] Our method, using a purely classical approach, is different from the aforementioned methods. It is theoretically simpler and more efficient in generating IR spectra with the MD simulations accelerated by GPUs. For example, it took us less than a day to calculate the spectra of Trpzip2 and the transmembrane domain peptide on using one Nvidia GTX1080 GPU. Most importantly, the main advantage of the method when compared to other methods is that our method allows for calculating transient IR spectra directly from nonequilibrium simulations. In other words, the laser-induced conformation changes, such as fibril dissociation and the IR spectra is directly obtained. Note that the direct equilibrium simulation of fibril dissociation is impractical with the current computer power, therefore the vibrational frequency maps cannot be used to obtain the transient conformations. Thus, our approach complements other methods. Nevertheless, the current theory which only considers C=O bond stretching vibration has limitations. First, although our method can reproduce the locations of the experimental absorption peaks well, the distribution widths, however, may be narrower than those found for the Tryptophan Zipper2 protein. Second, the anharmonicity effects were not completely considered because only implicit coupling between carbonyls was counted, while the contribution of the C-N bond absorption was neglected. We could observe the delocalization due to the natural coupling between C=O bond vibrations (Figure S1). Of note, the carbonyls were linked by the backbone and NH-CO hydrogen bonds in the MD simulation, resulting in a coupling effect between the C=O bond vibrations. Our computational algorithm can be further improved to take the full anharmonicity effect into consideration.

Applications of unpolarized laser all-atom MD simulations in phototherapies

Besides calculating IR spectra for a biomolecular system, our method may also be applied to find optimal conditions of phototherapies. The radiation absorption of biomolecular systems has been studied and widely applied for many health therapies such as infrared and visible light therapies.

The phototherapies have been medically used as the treatments for skin condition, cancer, sleep disorder, etc. [START_REF] Lee | UVB phototherapy and skin cancer risk: a review of the literature[END_REF][START_REF] Willis | A historical justification for and retrospective analysis of the systematic application of light therapy in Parkinson's disease[END_REF][START_REF] Hanford | Light Therapy and Alzheimer's Disease and Related Dementia: Past, Present, and Future[END_REF] Although the effectiveness of the phototherapies is evidently confirmed, the risks still exist. [START_REF] Lee | UVB phototherapy and skin cancer risk: a review of the literature[END_REF] Therefore, understanding the atomistic mechanism of the radiation absorption is vital for the therapy usage. It is still a grand challenge to fully reveal the atomistic picture of the phototherapies experimentally. Although computational approaches can describe biomolecular systems at the atomic level, a model/tool to mimic an unpolarized radiation in MD simulation is not available yet. Thus, the approach and implementation shown in this study will strongly support future studies to investigate the optimal conditions of phototherapies. Along this direction, we will focus on the development of the IR laser as a non-pharmaceutical approach for the treatment of Alzheimer's disease (AD). This is important because despite intensive studies in both academia and the pharmaceutical industry for many decades, there is still a lack of drugs, including small molecules, vaccines and antibodies, that could arrest or reverse the cognitive impairment of AD. A structural feature of amyloid fibrils in AD is the cross β-sheet stacking conformation stabilized by the hydrogen bonding network involved in the C=O bonds. In our previous studies, we have used the polarized IR laser to excite the amide I modes of the fibrils, resulting in the dissociation of the fibril structures [cite here our papers]. We find that the primary step in the dissociation process is due to the strong resonance between the fibril amide I vibrations and the tuned laser frequencies, and not just the deposited infrared thermal energy. Though a combination of experiments and simulations on disassembling of amyloid fibril of the 7-residue peptide GNNQQNY from yeast prion protein, we revealed the dissociation mechanism in unprecedented detail [Cite here our paper with Kawasaki]. Although the results are encouraging, the use of the polarized IR laser to target the amide I modes required high laser intensity in simulations compared to that used in experiments, such that the dissociation can be observed within reasonable simulation timescales. With the new unpolarized IR laser simulation method developed in this work, we will be able to directly target the C=O bonds. As a result, the required laser intensity should be lower, making it possible to directly compare simulations with experiments. This work is underway.

CONCLUSIONS

In conclusion, we have developed a comprehensive method for calculating the IR spectrum of the amide I modes of a single protein or amyloid fibril via all-atom MD simulations. The method involves mimicking an unpolarized laser that focuses the laser energy on the C=O bonds, which is equivalent to the effects of an IR polarized laser on an ensemble of identical proteins or a single molecule but rotating under a polarized laser. Technically, the method has been implemented into the AMBER simulation package and can be used with the fast GPU version of the pmemd.cuda module. We have tested the performance of the method with various single proteins and amyloid fibrils, including those with both β-sheet and α-helix structures. The simulated IR spectra are generally about 40 cm -1 blueshifted compared to experimental data. This is likely due to the inaccuracy of classical force fields which are not well-parametrized for high frequency bond stretching modes such as the C=O bond. We propose a simple correction by suitably rescaling of the fore constant of the C=O bond taken from the force fields. With this correction, we can produce not only the experimental spectra but also the difference between the peak positions of the α-helix and β-sheet structures. Moreover, the simulation method can also reveal small differences between the IR spectra of the parallel and anti-parallel β-sheet structures of proteins, consistent with previous high level theoretical calculations. We believe that the method will be useful for obtaining the amide I IR spectrum for proteins and amyloids, given its simplicity in theory, efficiency in computation, and accuracy in the result.

Our simulation method can also be combined with the IR pump-probe experimental technique\cite{PNAS 104, 14197, 2007} to track conformational changes in proteins and amyloids. In this experimental technique, an ultrashort IR laser pulse is used to excite the system, generating a non-equilibrium state, and another pulse is used to monitor the time evolution of the IR spectrum. In our simulation, the unpolarized laser excites the C=O bonds, resulting in conformational changes in the molecule. By calculating the IR spectrum for each conformation 20 during the simulation, we can obtain the time evolution of the spectra. Comparing the simulated spectra with experimental counterparts allows us to identify the conformations from the simulation trajectory whose IR spectra closely match those obtained from the experiment. We plan to develop this simulation/experimental approach to study the molecular mechanism of the IR laser-induced dissociation of amyloid oligomers and fibrils. This could contribute to the development of non-pharmaceutical approaches for treating Alzheimer's

ASSOCIATED CONTENT

Supporting information (SI) shows the details of the simulated systems (Table S1) and the spectra of dimer anti-parallel -sheet when the laser was employed to two residues and individual residue. (Figure S1)

  the amide I spectra with 11 values of dt (2, 1, 0.5, 0.4, 0.25, 0.2, 0.1, 0.08, 0.05, 0.02 and 0.01 fs). Figure4a, 4b and 4cshow the spectra of the different protein structures at dt = 1, 0.4, 0.2, 0.1 and 0.01 fs. The results indicate that when dt decreases, the peak of the spectrum first shifts

Figure 1 :

 1 Figure 1: Laser pulse mimicked by an external electric field and representative protein secondary structures studied in this work. (a) cartoon of the carbonyl (C=O) bond vibration (dotted blue box) when the unpolarized laser shown in the figure is applied, (b) representative protein secondary structures.

Figure 2 :

 2 Figure 2: The time dependence of the energy difference and the fluctuation of the C=O bond during the nonequilibrium MD simulations under an unpolarized laser for different protein structures: antiparallel β-sheet (a, b), parallel β-sheet (c, d) and α-helix (e, f). The energy difference of a system was calculated by dividing the difference between instantaneous energy and the energy at the beginning of the laser simulation by the number of targeted C=O bonds.The fluctuation was averaged overt all the C=O bonds. The laser's parameters are = 1.8 V/nm, = 3 ps, = 1 ps and different frequencies, . The result with laser frequencies of 1630, 1676, 1679, 1691 and 1720 cm -1 are shown in black, red, green, blue and orange colors, respectively. The data was averaged over ten independent MD trajectories.

Figure 3 :

 3 Figure 3: View of the external electric field with frequency, = 6000 cm -1 for different simulated time steps.

Figure 4 :

 4 Figure 4: The frequency dependence of carbonyl absorption energy in antiparallel β-sheet (a), parallel β-sheet (b) and α-helix (c) secondary structures at different simulated time steps. Panel (d) shows the dependence on simulation time step of the space frequency at which maximum absorption energy is observed.

Figure 5 :

 5 Figure 5: The effects of the laser parameters, E 0 and σ, on the absorption energy and structure of the α-helix model (PDB ID: 6N8C). The laser frequency is 1692 cm -1 which gives maximum absorption energy in the model. The result is similar for the β-sheet models (data not shown).

Figure 6 :

 6 Figure 6: The frequency dependence of the C=O bond absorption energy (a) and maximum fluctuation (b) in different protein secondary structures (PDB IDs: 2KIB, 1YJP, 6N8C). The data was obtained from the laser simulations with parameters, = 1.8 V/nm, = 1 ps and simulation time step dt = 0.1 fs.

Figure 7 :

 7 Figure 7: The frequency dependence of the C=O bond absorption energy antiparallel β-sheet (a) and α-helix (b) protein structures at different C=O force bond constants, . The = 570 is the original value of C=O force bond constant defined in the AMBER force field. The unit of the force bond constant is kcal/mol/Å 2 .

Figure 8 :

 8 Figure 8: IR spectra of the Trpzip2 peptide (PDB ID: 1LE1). The experimental data is from Ref.18. The calculated spectrum was normalized to the experimental value by matching the maximum absorption intensities of calculated and experimental spectra.

Figure 9 :

 9 Figure 9: IR spectra of the three different peptides encompassing the transmembrane domain of influenza A virus M2 H+ channel (PDB ID: 1MP6). The top, middle, and bottom panels are spectra from peptides that contain no isotopic label, a 13 C-16 O alanine residue at position 29 and 13 C-18 O glycine residue at position 34, respectively. The experimental data is from Ref.54. The calculated spectra were normalized to the experimental value by matching the maximum absorption intensities of calculated and experimental spectra.

ACKNOWLEDGEMENTS

This work was supported by the funds from the National Institutes of Health (NIH) NIH R01GM147673, NIH K25AG070277, and the National Science Foundation (NSF) 1955260. The authors thank the computing resources provided by the Center for Research Computing (CRC) at University of Pittsburgh.

AUTHOR INFORMATION Corresponding Author

Email: vhm3@pitt.edu and junmei.wang@pitt.edu

Author Contributions

V.H.M and J.W designed the project. V.H.M performed the simulations and analyzed data.

V.H.M, X.H, P.H.N, C.S, C.R, X.Q.X, and J.W discussed and wrote the paper.

Notes

The authors declare that they have no competing financial interests


Boris Guicheret

Laurent Vanoye

Camila Rivera-Cárcamo

Claude De Bellefon

Philippe Serp

Régis Philippe
email: regis.philippe@lgpc.cpe.fr

Alain Favre-Réguillon
email: alain.favrereguillon@lecnam.net

Solvent-free Hydrogenation of Squalene using Parts per Million Levels of Palladium Supported on Carbon Nanotubes: Shift from Batch Reactor to Continuous-Flow System

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

Introduction

The most challenging problem currently facing the pharmaceutical and specialty chemical industries is the transformation of existing batch catalytic processes into continuous flow processes. Such transformation promises many benefits, such as increased safety, simpler scale-up, access to novel processing windows, and higher throughput. This transition could only be possible if highly active and, as importantly, stable catalysts are developed. In addition, from a green chemistry perspective, solvent-free chemical transformations are more interesting than reactions carried out with any solvents. [1] However, from a process chemistry perspective, solvent-free reactions are linked to safety issues (especially for highly exothermic reactions), mass transfer limitation (for viscous products), and inhibition (particularly for bio-based substrates). [2] Therefore, the development of environmentally friendly solvent-free reactions continues to be of great interest.

The reduction of squalene (SQE) into squalane (SQA) is a representative illustration of this particular challenge (Figure 1). SQE is a linear triterpene present in the cells of all living organisms, which is extensively used in pharmaceuticals, food supplements, and cosmetics due to its multiple functions. [3] The global SQE market size was higher than USD 110 million in 2016 and is expected to witness significant growth over the following years. [3] In order to increase its stability (SQE oxidizes quickly due to its unsaturated character), SQE could be fully hydrogenated to SQA (Figure 1). Besides applications in cosmetics and pharmaceuticals, SQA, because of its high boiling point, inertness, low viscosity, and high gas solubility is an interesting apolar solvent used as the reaction medium for catalytic tests [4] or the synthesis of quantum dots. [5] Figure 1 Hydrogenation of squalene (SQE) into squalene (SQA)

The highly exothermic hydrogenation of SQE is traditionally carried out in stirred tank reactor using high loading of metal-based catalyst (Ni/diatomaceous earth, [6] copper chromite, [7] and Pd [1c, 1d, 8] ) under harsh conditions (Table 1 for Pd catalysts). Under those conditions, fast and irreversible catalyst deactivation and metal leaching occur, requiring product purification and catalyst refining and remanufacturing. [9] As such, purification by distillation could be necessary to obtain high-quality SQA, which increases overall production costs. Furthermore, safety issues arise due to the exothermicity of the reduction (HR = -765 KJ.mol -1 ). Without solvent, an adiabatic temperature rises higher than 860 K to reduce SQE could be calculated. Thus, intentionally degraded operation conditions for safety reasons are generally implemented at an industrial scale, i.e., a two-step, one-pot process. The first step is performed at low H2 pressure to initiate the reduction without the risk of a thermal run-away, followed by a subsequent step at higher H2 pressure and temperature to obtain the total reduction (Table 1, entries 1 to 3). At the lab scale, solvent-free SQE could be fully hydrogenated using new Pdbased catalysts, but prolonged reaction time and harsh conditions are still required (Table 1, entries 3 and 4). Clearly, there is a need for more active and stable hydrogenation catalysts to convert SQE into SQA under milder conditions, especially now that bio-based SQE is poised to entirely replace SQE of animal origin. 3 Pd/clay 1 300 °C [b] -10 bar [c] This work [a] Not provided. The reaction was stopped when H2 uptake ceased. [b] Temperature was gradually increased from 25 °C at a rate of 5 °C.min -1 [c] H2 pressure at 25°C

Recently, we have prepared new Pd-supported catalysts on carbon nanotubes (CNT), bearing both Pd nanoparticles (PdNP) and Pd single atoms (PdSA). [10] We have demonstrated that the control of the ratio between PdSA and PdNP on carbon supports enables, thanks to cooperative catalysis, [11] the development of very active hydrogenation catalysts integrating the ultra-rational use of precious metals. [12] In particular, we have reported that Pd/CNT catalysts show a high activity compared to catalysts already described for hydrogenation of terpenes, i.e., myrcene and SQE in solution (0.33 M in heptane). [12b] These catalysts should meet the need to increase catalytic activity for a highly demanding reaction like solvent-free hydrogenation of SQE.

Herein, we compared the Pd/CNT catalyst activity with commercial catalysts, and we evaluated the Pd/CNT catalyst stability in continuous flow reactors after its immobilization on a metallic foam. Then, an intensified commercial pilot-scale continuous flow reactor was implemented to demonstrate its efficiency and robustness and provide a safe, compact, intensified, and scalable production process developed for the solvent-free hydrogenation of SQE.

Results and Discussion

Solvent-free hydrogenation of SQE with commercial catalysts

The solvent-free reduction of SQE was first studied using commercial catalysts as a benchmark. The reaction progress could not be easily followed by in-situ or ex-situ analysis. SQE has 6 double bonds with similar reactivity (all tri-substituted double bonds), and partial reduction leads to the formation of many isomers and diastereoisomers (there are 4 prochiral carbon in SQE). Nevertheless, reactant (SQE) and product (SQA) could be identified and quantified from the bulk of intermediates using GC/FID using tetradecane as an internal standard (Figure 2). Furthermore, monitoring H2 consumption as the function of time (

Figure 3), associated with off-line GC/FID analysis, could give highly reproducible monitoring of the catalyst activity. By using commercial 5wt % Pd/-Al2O3, under typical conditions (120 °C, H2 20 bar, 2 mol% Pd), the solvent-free catalytic reduction of SQE was achieved after 8.6 h (Figure 3 a). The time required to achieve the total transformation of SQE into SQA is comparable to published data (Table 1). For SQE hydrogenation in heptane (0.33 M), complete saturation was obtained in less than 0.3 h for equivalent experimental conditions (Figures S5 andS6). The strong adsorption of SQE could explain this concentration effect and partially hydrogenated products on the catalyst surface, resulting kinetically in a term of inhibition by the reactant on the Langmuir-Hinshelwood-Hoogen-Watson (LHHW) rate law. [1d, 13] The reduction of solvent-free SQE was also studied with 2 mol % Ni using a commercial 65wt % Ni/SiO2-Al2O3 catalyst under the same conditions (Figure 3 b). After 30 min at 120 °C, no H2 consumption could be detected, and thus the temperature of the reactor was heated up to 150 °C. Signal (a.u.)

Then a slow reduction of the solvent-free SQE could be observed, and total reduction could be obtained after more than 5 h at 150 °C. Solvent-free hydrogenation of SQE with Pd/CNT.

Our groups recently reported that Pd supported on mesoporous nano-carbons catalysts showed high activity for the hydrogenation of terpenes in heptane (-myrcene, SQE) [10,12] and demonstrated that PdSA and PdNP were simultaneously present on the carbon supports. Unique cooperative catalysis operates between these species, allowing one to reach very high activity for alkene hydrogenation through an H-spillover assisted mechanism involving the support functionalized surface. [12a, 14] On CNT support, the optimum PdSA/PdNP ratio was found at 10/1 for terpene hydrogenation in solution. [12] The main characteristics of the Pd/CNT catalyst used for SQE hydrogenation are summarized in the experimental section (other characterization can be found in ref. [12a] ), and Figure 4 shows representative STEM-HAADF micrographs (PdSA are circled in yellow for better visualization). It has been previously shown that the activities of such Pd/CNT catalysts are very high. [12] For safety reasons, the amount of Pd was thus decreased to 0.03 mol % (300 ppm). The solvent-free reduction of SQE was first studied under standard conditions, i.e., 120 °C, 20 bar H2 (Figure 5). Under those relatively mild conditions, the total and selective reduction of SQE into SQA was achieved in less than 1.4 h with 300 ppm of Pd (Table 1, entry 6). Despite using an efficient Ranque-Hilsch vortex device for cooling the relatively small stirred tank reactor (200 mL), the liquid phase temperature reaches 130°C within the first min of reaction. Such heat release obtained for 300 ppm of Pd prevents the use of higher loading of catalyst in a batch reactor. However, this result clearly demonstrates that Pd supported on CNT that combines PdSA and PdNP can reach a very high activity for SQE hydrogenation under solvent-free conditions. This catalyst is much more efficient than the catalysts already described in the literature (Table 1) for the solvent-free reduction of SQE. The implementation of this catalyst in a continuous flow reactor was performed to evaluate its stability. Furthermore, the use of an intensified reactor with efficient heat transfer will enable us to study new process windows for this reaction [15] and to propose a more intense/efficient reactor for heat transfer and thus be able to propose an intensified synthesis with new operating conditions

Assessment of catalyst stability

The catalyst stability could be evaluated by several consecutive batch reactions under the same conditions. In each cycle, after completion of the reduction, the catalyst should be collected by filtration, washed, dried, and reused in a subsequent cycle. Under those conditions (Figure 5), a maximum Turn Over Number (TON) of 20'000 could be achieved for such a cycle.

Otherwise, implementing this catalyst in a continuous flow reactor is an outstanding approach to demonstrate heterogeneous catalysts' stability with access to very high TON. However, implementing Pd/CNT catalyst in a continuous flow microreactor is not straightforward since Pd/CNT is a nanotextured micro powder (particle size < 10 m). Heterogeneously catalyzed reactions with powder catalysts typically used in stirred tank reactors (particle size between 10-100 m) have been successfully used in microreactors. [16] The catalysts could be confined in micro-packed bed systems [17] , but limitations such as high-pressure drops and possible temperature gradients are likely. This is especially relevant with the viscous SQE/SQA charge (36.1 and 2.2 mPa.s -1 for SQA at 20 °C and 120 °C compared to 0.39 mPa.s -1 for heptane at 20°C). [18] New strategies, adapted explicitly to micro-or milli- reactors, have been developed, such as catalytic static mixers (CSM), [19] and metallic open cell foams (OCF). [20] These scaffolds, coated with an active porous catalyst layer, can be inserted into continuous flow reactors. CSM should be specifically designed to fit in the microreactors and could be produced by additive manufacturing. [19b, 19c, 21] Another approach would be to use high voidage structured materials commercially available at low cost, i.e., OCF. With their excellent thermal performance and large external specific surface area (up to 10,000 m 2 /m 3 ), OCF have demonstrated their interest as heat exchangers, internal mixers, and catalyst supports. [20,22] Fecralloy metallic foams (Hollomet, Foamet 400µm, void fraction 92 %, and geometrical specific surface area of 7000 m 2 /m 3 ) could be precisely adapted as milli-reactors internals from a foam sheet by electro-machining discharge. This offers a good cutting accuracy (±5 m) for conductive substrates, ensuring adaptability to various reactor geometries and intimate contact with channel walls, avoiding by-pass flows. OCF millimeter-scale blocks were then dip-coated with an aqueous suspension of the Pd/CNT catalyst (100 g/L) according to a reported procedure. [23] A final heat treatment at 300 °C under a reductive N2-H2 mixture provides the final supported catalyst. The coated OCF were controlled by optical microscopy (Figure S 7) and elemental analysis. A mean catalyst content of 16.5  5.0 g Pd/CNT catalyst per L of foam was obtained. The coating process did not significantly modify the particle size distribution (2 ± 0.5 nm) and the PdSA/PdNP ratio (5.6) of the Pd/CNT catalyst (Figure S 8), pointing to its high thermal stability.

The coated OCF were inserted into a homemade millimeter-scale tubular reactor (4 tubes of 100 mm length and 4.4 mm internal diameter) to provide an effective reactor volume (VR) of 5.8 mL with a total catalyst amount of 0.1 g. These tubes are positioned in series in a thermoregulated oven (Figure 6), and thermocouples distributed along the tube allow for measuring the axial temperature profile (See Supporting Information for details). The continuous experiments were performed at 180 °C and 30 bar of H2 with solvent-free SQE (2.03 M) to obtain a total reduction of SQE. An H2/SQE molar ratio higher than 14 was used to ensure a large excess of H2. First, the influence of contact time, defined as the ratio between the effective reactor volume (VR) and the liquid flow (QL), on the yield in SQA and conversion of SQE was studied in this configuration. The results are shown in Figure 7. For a contact time below 6 min at 180 °C and 30 bar H2, total conversion of SQE was obtained, and the yield in SQA was up to 70 % (Figure 7). These results further emphasized the synergy between the high activity of the Pd/CNT catalyst and the mixing effect of the metallic foam on H2 mass transfer. No increase in the liquid phase temperature inside the reactor could be noticed. In order to increase the yield in SQA, the contact time was increased by decreasing the liquid flow (QL). For a contact time higher than 14 min, at 180 °C and 30 bar of H2, the yield of SQA was higher than 98%.

Figure 8 depicts 2 distinct curves to illustrate the stability along time on stream in this continuous reactor. The first one (squares) reflects the entire catalyst work through the evolution of TON in the continuous hydrogenation of SQE. It is noticeable that this curve is obtained under various experimental conditions, i.e., temperature, gas and liquid flow, and concentration of SQE, explaining the variations in slope. Under those various conditions, a cumulative turnover number (TON) up to 300'000 molH2/molPd could be determined (Figure 8). It should be mentioned that such a value is largely minimized compared to the real value since only the amount of SQA produced was considered to calculate the TON (eq. 1). Intermediates, i.e., partially hydrogenated SQE, should theoretically be taken into account for the determination of the TON. However, since they could not be precisely quantified (Figure 2), TON was only based on the SQA produced (eq. 1).

𝑇𝑂𝑁 =

6×mol of SQA produced mol of Pd (eq. 1)

The second curve (triangles) is dedicated to periodic control experiments of the space-time yield using rigorously identical experimental conditions (i.e., SQE 0.3 M in heptane, contact time 2.9 min, 120 °C, and 30 bar of H2). These activity control experiments were performed regularly throughout the different conditions experienced by this catalytic object in the continuous reactor. Under those conditions, a space time yield (STY-90%), defined by the quantity of SQA produced in the reactor per second and per mol of palladium, of 36  5 molSQA.molPd -1 .s -1 was obtained. No decrease of the STY-90% as the function of time on stream could be noticed for a time on stream higher than 18h (Figure 8). The STEM-HAADF analyses carried out on the spent catalyst do not show any significant modification in the size of the palladium particles (2.4 ± 0.1 nm) and the PdSA/PdNP ratio (6.1) (Figure S9), pointing to its high thermochemical stability.

Scale-up and transposition to a commercial intensified pilot reactor

With highly efficient catalysts, demanding batch hydrogenation protocols pose safety hazards, particularly when large-scale reactors are required, because of reduced heat transfer capacities compared to lab-scale. Furthermore, the resulting temperature excursions may affect the selectivity of the transformations and cause the deactivation of the catalysts. This is particularly significant in the case of solvent-free SQE hydrogenation, where an adiabatic temperature rise of 863 K could be estimated. These safety concerns associated with mass transfer-limited reaction kinetics hindered the further development of highly efficient hydrogenation catalysts in physically limited reactors.

Contrary to stirred tank batch reactors, heat exchange in structured micro-/milli-reactors is sufficiently intense that even highly exothermic reactions can be conducted under nearly isothermal conditions. [24] Furthermore, due to the large interfacial areas, very efficient hydrogen mass transfer is achieved. [25] Finally, for a given production capacity, continuous processes are inherently safer because of drastically reduced reactor inventories in opposition to stirred tank batch reactors. For these reasons, continuous-flow hydrogenation protocols employing metal catalysts have increased significantly over the past few years. [16a, 26] In order to demonstrate the robustness and the scalability of the results obtained in the homemade laboratory reactor, a transposition of the solvent-free SQE hydrogenation to a commercial intensified reactor was studied. The Miprowa-Lab reactor from Ehrfeld Mikrotechnik [27] (Figure 9a) was chosen for its similar millimeter-range hydraulic diameter (12 x 1.5 mm 2 ) and its possible and easy operation with coated foam internals (Figure 9b). This reactor presents a practical reacting volume of 43.2 mL (divided into 8 channels in series), corresponding to a potential scale-up factor of 7.5 compared to our homemade tubular reactor (5.8 mL). Again, foam pieces were precisely cut in the same foam block to be mounted inside the 8 channels after coating with the same Pd/CNT catalyst (Figure 9c). A catalyst content of 16.5  4.0 g of Pd/CNT catalyst per L of foam was obtained, close to the one obtained for the tubular reactor. Table 2 compares the two reactors for two different contact times, i.e., 14.4 min and 9.6 min, under the same experimental conditions (180 °C, 30 bar H2, solvent-free SQE). For a given contact time, the SQA yields in different reactors are close, demonstrating the scalability of the coated foam contacting mode to a commercial reactor. A slight difference could be noticed, reflecting both the slight difference in catalyst content and possible variations in heat and mass transfer capacities due to the changes in the geometrical configuration of the reacting volume (channel section, channel length, aspect ratio, etc…) and wall heat exchanging system.

Transposition of SQE reduction at two similar operating points, from the homemade tubular reactor to the commercial intensified reactor, results in higher production capacity (Table 2). From lab-scale tubular reactor to Miprowa-Lab reactor, the production capacity of high purity SQA (>98 %) was increased by a factor of 7.5, i.e., 0.485 kgSQA/day to 3.630 kgSQA/day for tubular reactor and Miprowa-Lab reactor, respectively. Under those conditions, even if all Pd present in the microreactor was leached in 24h, the amount of Pd present in SQA will be less than 50 ppb, far below the parenteral concentration limit for Pd in drug products, i.e., 1 ppm (1 μg/g) due to the toxicity of Pd (ICH classification 2b). [START_REF]ICH Harmonised Guideline, Guideline for Elemental Impurities, Q3D(R1) 22 march 2019, International Council For Harmonisation of Technical Requirements For Pharmaceuticals For Human Use[END_REF] It should be noticed that the absence of appreciable deactivation during the 18h time on stream study (Figure 8) indicates that the situation concerning Pd leaching is far from this most pessimistic scenario. SQA yield [b] [%] STY [molSQA/molPd/min] SQA Production capacity [c] [kgSQA/day] [a] Reaction conditions: solvent-free SQE (2.03 M) with tetradecane as internal standard (0.1 M); 180 °C; H2 30 bar; Pd/CNT coated on Fecralloy metallic foams (Foamet 400µm). Catalyst content of 17.2 g and 16.5 g of Pd/CNT per L of foam for tubular and Miprowa-Lab reactor, respectively. QG = 400 NmL/min and 3000 NmL/min for tubular and Miprowa-Lab reactors, respectively. [b] SQA yield was determined by GC/FID using internal calibration. [c] SQA purity was higher than 98% without could be used without further purification

In case of coated foam availability, further seamless scale-up would be quickly and easily ensured (by the technology provider) up to reactor units able to process several m 3 /h. [27] Other internals could also be coated, such as the pristine static mixing elements of Ehrfeld microreactors or CSM. [19b, 19c] Before such demonstrations, identical transfer features need to be demonstrated, and the Pd/CNT catalyst coating has to be adapted. Furthermore, except for the first double bond hydrogenation of SQE, the reactor operates roughly in a chemical regime. Further process intensification can be expected thanks to this equipment with higher H2 pressure and a slightly higher temperature operating window, and it still needs to be investigated.

Conclusion

The catalytic activity of palladium catalyst with PdSA and PdNP supported on CNT was tested on solventfree squalene's hydrogenation. The catalytic activity of Pd/CNT showed a tremendous improvement compared to previously described and/or commercially available catalysts. The total reduction of solvent-free SQE into SQA was effective at 120 °C, 30 bar H2 using 300 ppm of Pd in less than 1.4 h in a lab-scale batch stirred tank reactor. The catalyst stability was assessed in a continuous flow reactor, and a TON higher than 300'000 without any detectable loss in activity was demonstrated. A successful batch to a commercial continuous milli-structured pilot-scale reactor was then validated. These experiments assess the operability of such a catalyst in flow and its stability. Neat squalane ( 98%) could be obtained continuously at 180 °C and 30 bar H2 with a contact time below 15 min. Under those conditions, a space-time yield of 43.7 molSQA.molPd -1 .min -1 and a high purity squalane production capacity of 3.6 kg per day was obtained using a commercial intensified reactor with an effective reactor volume (VR) of 43.2 mL.

Experimental Section

Squalene and solvents SQE (98%) was purchased from Aldrich and used as received. SQE was stored under refrigeration (2-8°C) and protected from light. Heptane, EtOH, Triton X-100, and dextrin were purchased from Aldrich or Acros and used as received.

Catalysts

Pd 0 /Al2O3 (Alfa Aesar, Product number 11713) and Ni/SiO2-Al2O3 (Alfa Aesar, Product number 31276) were used directly without further treatment. Pd/CNT was prepared by wet impregnation method according to P. Serp et al.. [10] The main characteristics of these catalysts are summarized in Table 3, and Figure 4 shows representative STEM-HAADF. Particle size [a] [nm]

Pd dispersion [b] [%]

PdSA/PdNP ratio [c] Pd/Al2O3 [d] 5.00 140 n. -Pd/CNT [e] 2.04 232 1.5 ± 0.5 67 3.1 [a] from TEM [b] Metal dispersion was evaluated from a universal mathematical relation between the mean relative size of metallic crystallites and their dispersion, as described elsewhere ref. [START_REF] Borodziński | [END_REF] . [c] The SA/NP ratio is a number ratio measured from the STEM-HAADF analyses [d] see ref. [30] [e] see ref. [10, 12a]

Characterization of the catalysts

The palladium content in the products was measured by inductively coupled plasma optical emission spectroscopy (ICP-OES). TEM and HR-TEM analyses were performed at the "Centre de microcaracterisation Raimond Castaing, UMS 3623, Toulouse" by using a JEOL JEM 2100F equipped with a field emission gun (FEG) operating at 200 kV with a point resolution of 2.3 Å and a JEOL JEM-ARM200F Cold FEG operating at 200 kV with a point resolution of > 1.9 Å. The particle size distribution was determined through a manual measurement of enlarged micrographs from different areas of the TEM grid (at least 300 particles). The size distribution reported as Pd atomic percentage (based on total atoms) in each size range was obtained by calculating the number of atoms (N) in each nanoparticle, assuming a spherical shape using an equation described in ref. [31] . The SA/NP ratio (a number ratio) was measured from the STEM-HAADF analyses of at least 500 elements. Optical microscopy was performed using optical microscope KEYENCE VHX 6000.

General procedure for the hydrogenation of SQE in a stirred tank reactor

The reaction was performed in a 200 mL thermo-regulated stainless steel autoclave from Top Industrie. A stainless steel insert with four baffles was mounted on the reactor walls, and the reaction medium was stirred by a gas-inducing Rushton impeller having six straight blades. Temperature control is provided by both an electric shell (for heating) and a Ranque-Hilsch vortex device (for cooling). The reactor is operated in batch mode for the liquid and solid phases and semi-batch mode for the gas phase. Hydrogen was constantly fed via a pressure regulator and supplied by a calibrated reserve. The pressure and temperature monitoring of this reserve provide the hydrogen consumption rate. The liquid phase temperature, monitored using J type thermocouples, and the time-course of hydrogen consumption, monitored with a pressure transducer, were acquired and recorded online using Labview® at 1 Hz.

The catalyst and 80 mL of solvent-free SQE, with tetradecane (0.1 M) as an internal standard, were placed in the reactor. The reactor was purged three times with N2 at room temperature, and the reactor was heated up to the desired temperature under stirring. When the desired temperature was reached, stirring was stopped. The N2 atmosphere was replaced by the H2, and the experiment was restarted by switching on the agitator, determining the t0 of the reaction.

At the end of each experiment, the reaction mixture was analyzed after dilution with heptane using an Agilent 6890 GC-FID and a Shimadzu 2010 GC-MS system, both equipped with DB-17 columns. Unlike other terpenes, the reactivity of the double bonds in SQE is quasi-identical (all tri-substituted double bonds). As such, the identification of different isomers as the function of time by classical GC and GC/MS methods is impractical. Nevertheless, reactant (SQE) and product (SQA) could be identified and quantified from the bulk of intermediates using GC/FID using tetradecane as the internal standard (Figure 2). In-line H2 consumption associated with off-line GC analysis could give highly reproducible catalyst activity monitoring.

Dip-coating procedure

Preparation of the suspension: In the tank of a planetary mill, 20 mL of water was added with 1.1 mL of Triton X-100 and stirred at 300 rpm for 3 min. 0.8 g of dextrin was added, and the mixture was ballmilled for 5 min at 300 rpm. 5g of Pd/CNT was added, and the mixture was again ball-milled for 45 min at 500 rpm. The final suspension of the catalyst was collected with water to reach a final concentration of 100 g/L. Activation of the supports: Fecralloy metallic foams (Hollomet, Foamet 400µm) were cut into cylinder-shaped blocks (length-diameter: 50 mm-4.3mm) or into sheets (length, width, thickness: 150 mm, 12 mm, 1.5 mm). The metallic foams were carefully washed in acetone in a sonication bath, then dried and heated to 500 °C for 4 hours under 500 NmL/min airflow and then cooled to room temperature. Dip-coating of the support: Activated substrates were dipped entirely in the catalyst suspension. Excess catalyst suspension was removed by airflow. The coated foams were dried under air at 80 °C overnight, then heated to 300 °C for 2 hours under nitrogen/hydrogen flow (80/20) (total flowrate of 500 NmL/min) and kept for 2 hours at 300 °C. The oven was cooled down to room temperature under nitrogen flow (500 mL/min). The amount of catalyst coated on the metallic foams was determined by weighing. Analysis of the Pd/CNT catalyst after coating shows that the coating process did not significantly modify the particle size distribution (2 ± 0.5 nm) and the PdSA/PdNP ratio (5.6) (Figure S 8).

General procedure for the hydrogenation of SQE in flow in a milli-reactor Coated foam pieces were inserted in a homemade milli-reactor (4 tubes in series of 100 mm length and 4.4 mm internal diameter) or in a Miprowa-Lab reactor (8 channels in series of 300 mm length with a width of 12 mm and a height of 1.5 mm). The temperature was controlled by a thermoregulated oven or a heated circulating bath (Figure 6 and Supporting Information). The temperature inside the reactor was monitored using J-type thermocouples.

In a 1000 mL volumetric flask, solvent-free SQE with tetradecane as internal standard (0.1M), was prepared. At the reactor inlet, solvent-free SQE (ensured by an HPLC pump Shimadzu LC 20 AD) and gas flow (regulated by a mass flow controller Bronkhorst Elflow Prestige) meet at a T junction. An Equilibar back pressure regulator precisely regulated outlet pressure. G-L separation and liquid collection are ensured by gravity. When a stable temperature was obtained, periodical analysis of the liquid phase by GC-FID and GC-MS after dilution with heptane was used to calculate the conversion of SQE and the yield of SQA.

Table of Contents

Solvent-free total hydrogenation of highly unsaturated squalene into squalane is achieved in batch in less than 2h using Pd-supported catalyst on carbon nanotubes bearing both Pd nanoparticles and Pd single atoms under mild conditions (120°C, 20 bar, 300 ppm Pd). Implementation of this catalyst in a commercial continuous flow pilot was successfully achieved and a TON higher than 300'000 was demonstrated.

H 2 solvent-free

Figure 2 .

 2 Figure 2. Comparison of the GC/FID chromatograms as the function of reaction progress. a) initial composition of reaction mixture; b) reaction mixture at 85% conversion, 1 h; c) final composition of the reaction mixture after 9 h (Conversion >98%). Reaction conditions: H2 20 bar; 120 °C; liquid phase 80 mL of solvent-free SQE (2.03 M) with tetradecane as internal standard (0.1 M); Pd 0.08 mol % (280 mg of 5 wt % Pd/-Al2O3), 120 °C. Heptane was used to dilute the reaction mixture for GC analysis. (See supporting information in Figures S1 to S4 for full-scale chromatograms).

Figure 3 .

 3 Figure 3. H2 consumption and temperature of the liquid phase as the function of time. Reaction conditions: H2 20 bar; 120 °C; liquid phase 80 mL of solvent-free SQE (2.03 M) with tetradecane as internal standard (0.1 M); a) Pd 0.08 mol % (280 mg of 5 wt % Pd/-Al2O3), 120 °C; b) Ni 2 mol % (300 mg of 65 wt % Ni/SiO2-Al2O3), 120 °C for 0.5 h then 150 °C.

Figure 4 .

 4 Figure 4. STEM-HAADF micrographs and particle size distribution are based on total particles and total atoms for Pd/CNT catalyst.

Figure 5 .

 5 Figure 5. H2 consumption and temperature as the function of time. Reaction conditions: H2 20 bar; 120 °C; liquid phase 80 mL of solvent-free SQE (2.03 M) with tetradecane as internal standard (0.1 M); Pd 0.03 mol % -300 ppm (270 mg of 2.04 w% Pd/CNT).

Figure 6 .

 6 Figure 6. Experimental setup for the assessment of Pd/CNT catalytic stability.

Figure 7 .

 7 Figure 7. Composition at the outlet of the reactor as the function of contact time. Reaction conditions: effective reactor volume (VR) 5.8 mL; 95.7 mg of 2.04 wt% of Pd/CNT coated on Fecralloy OCF (Foamet 400 µm); solvent-free SQE (2.03 M) with tetradecane as internal standard (0.1 M); 180 °C; H2 30 bar; (QG/QL)0 >30. SQA and SQE concentrations were determined by GC/FID using internal calibration. The concentration of the different intermediates was calculated by the difference in the concentration of SQA and SQE (see text and Figure2). Lines are guides to the eyes.

Figure 9 .

 9 Figure 9. a) Ehrfeld Mikrotechnik Miprowa-Lab reactor; b) head of the reacting zone showing the foam sheets internals; c) coated foam sheets with Pd/CNT catalyst.

Table 1 .

 1 Experimental conditions for solvent-free hydrogenation of SQE by Pd catalysts.

	Entry	Catalysts	Pd [mol%]	Reaction conditions [temp. -P H2]	Time [h]	Ref.
	1	Pd/C	0.1	120 °C -50 bar then 160 °C -150 bar	[a]	[8a]
	2	Pd/C	0.1	25° C -65 bar then 60 °C -65 bar	47	[8b]

Table 2 .

 2 Comparison of a homemade tubular reactor (VR = 5.8 mL) and commercial intensified reactor (VR = 43.2 mL) at two operating points[a]

	Effect	Contact time	QL
	[mL]	[min]	[mL/min]

Table 3 :

 3

	Catalysts' main characteristics
	Catalyst	Pd content	BET surf. area
		[wt%]	[m 2 .g -1 ]

Acknowledgments

Part of this work was performed within the framework of the "DEEPER" project funded by the region Auvergne-Rhône-Alpes (contract number 15 021131 01 -CNR006) through the 20 th FUI call. Part of this work was also supported by the Agence Nationale de la Recherche (ANR project ANR-19-CE07-0030, COMET), which is gratefully acknowledged.

Frédéric Bornette and Fabrice Campoli are acknowledged for contributing to the reactor design, fabrication, and instrumentation.

Ehrfeld Mikrotechnik was greatly thanked for the loan of the Miprowa-Lab reactor (www.ehrfeld.com/home.html).

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Supporting Information

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.