Petteri A Vainikka

Sebastian Thallmair

Paulo C T Souza

Siewert J Marrink
email: s.j.marrink@rug.nl

Martini 3 coarse-grained model for type III deep eutectic solvents: Thermodynamic, structural, and extraction properties

Deep eutectic solvents (DESs) are a more environmentally friendly, cost-effective, and recyclable alternative for ionic liquids. Since the number of possible deep eutectic solvents is very large, there are needs for effective methods to predict the physicochemical nature of possible new deep eutectic solvents that are not met by the currently available models. Here, we have built coarse-grained models for a few well-known and actively studied deep eutectic solvents using the recently published Martini 3 force field.

Molecular dynamics simulations demonstrate that our models predict the properties of these particular DESs with an acceptable accuracy, and are capable of capturing known liquid-liquid extraction processes as well as morphological shape changes of surfactant aggregates. Our coarse-grained approach is novel in the study of DESs, opening new possibilities for rapid screening of new DESs and their properties.

Introduction

Deep eutectic solvent (DES) is a name for a class of solvents which are typically formed as a result of mixing Brønsted or Lewis acids with hydrogen bond acceptors (HBA). [START_REF] Abbott | physical and electrolytic properties of eutectic mixtures based on CrCl3•6H2O and urea[END_REF][START_REF] Smith | Deep Eutectic Solvents (DESs) and Their Applications[END_REF] Within a correct molar ratio, which varies depending on the thermodynamic properties of the pure compounds forming the DES, they form an eutectic mixture, manifested by a significant depression in the freezing point, allowing them to remain liquid at room temperature. [START_REF] Martins | Insights into the Nature of Eutectic and Deep Eutectic Mixtures[END_REF] Although they are primarily synthesized for laboratory and industrial use, they are also present in nature: Recent work suggests that so called 'natural DESs' (NADESs) could form a separate liquid phase inside mammalian and plant cells, contributing to both the cellular metabolism, and the survivability of the cell. [START_REF] Choi | Are Natural Deep Eutectic Solvents the Missing Link in Understanding Cellular Metabolism and Physiology?[END_REF] Our work focuses primarily on type III DESs, which are formed by mixing a quaternary ammonium salt with a hydrogen bond donor. [START_REF] Smith | Deep Eutectic Solvents (DESs) and Their Applications[END_REF]5 Care should be taken when considering whether or not a solvent is indeed a 'deep eutectic' solvent: The simple presence of a eutectic point, which is ubiquitous in mixtures, or the ability to form hydrogen bonds is not enough. [START_REF] Martins | Insights into the Nature of Eutectic and Deep Eutectic Mixtures[END_REF]6 The work of Martins et al. [START_REF] Martins | Insights into the Nature of Eutectic and Deep Eutectic Mixtures[END_REF] is recommended for further consideration of the topic.

Publications related to using eutectic mixtures as substrates in enzymatic catalysis can be found in mid-1990s, 7,8 but the wide range of their potential uses as solvents was not realized until the early 2000s. In 2003 Andrew Abbott and colleagues reported the synthesis of a DES by mixing choline chloride and urea. 9 This is, to the extent of the authors knowledge, the first use of the term 'deep eutectic', coined in order to differentiate between ionic liquids (ILs) and DESs. In the following year Abbott et al. demonstrated the synthesis of multiple DESs using quaternary ammonium salts and carboxylic acids. 10 This work was among the first to highlight multiple, potentially important and lucrative properties of DESs: They share many of the coveted properties of ILs, such as wide liquid ranges, very low vapour pressures, non-flammability, and applicability in various processes relevant to both academia and industry. 6,[11][12][13] Contrary to ILs, preparation of DESs has optimal atom economy, 13 they require comparatively little purification as the purity of any DES is directly ordained by the purity of the compounds from which it was generated from, [START_REF] Martins | Insights into the Nature of Eutectic and Deep Eutectic Mixtures[END_REF]14 and most materials used in their formulation are cheap and readily available. [14][15][16] Additionally, by carefully selecting the precursors, DESs can be made environmentally friendly 17 and/or easily recyclable. 18 Although the total number of possible DESs is unknowable, and will likely to remain so until the community has come to an agreement on how to classify a solvent as a DES, current works have demonstrated that the growing amount of existing DESs is turning comprehensive experimental mapping of their properties and behaviour less and less feasible. 19,20 This has prompted the use of theoretical and computational methods in studying the properties and predicting new possible DESs. These methods include, but are not limited to, molecular dynamics simulations, [START_REF] Kaur | Anatomy of Microscopic Structure of Ethaline Deep Eutectic Solvent Decoded through Molecular Dynamics Simulations[END_REF][START_REF] Kumari | Influence of Hydration on the Structure of Reline Deep Eutectic Solvent: A Molecular Dynamics Study[END_REF][START_REF] Atilhan | Molecular dynamics simulations of mixed deep eutectic solvents and their interaction with nanomaterials[END_REF][START_REF] Triolo | Liquid structure of a choline chloride-water natural deep eutectic solvent: A molecular dynamics characterization[END_REF][START_REF] Salehi | Computing solubility parameters of deep eutectic solvents from Molecular Dynamics simulations[END_REF][START_REF] Kaur | Microstructure of Deep Eutectic Solvents: Current Understanding and Challenges[END_REF] classical thermodynamic models such as group contribution methods, 19,[START_REF] Shahbaz | Prediction of deep eutectic solvents densities at different temperatures[END_REF][START_REF] Adeyemi | Physicochemical properties of alkanolamine-choline chloride deep eutectic solvents: Measurements, group contribution and artificial intelligence prediction techniques[END_REF] and quantum chemical calculations. [START_REF] Bezold | Assessing solute partitioning in deep eutectic solvent-based biphasic systems using the predictive thermodynamic model COSMO-RS[END_REF][START_REF] Liu | Screening Deep Eutectic Solvents for CO2 Capture With COSMO-RS[END_REF][START_REF] Abranches | Using COSMO-RS to design choline chloride pharmaceutical eutectic solvents[END_REF][START_REF] Gutiérrez | Theoretical Study on Deep Eutectic Solvents as Vehicles for the Delivery of Anesthetics[END_REF] While some group contribution methods (COSMO-RS, UNIFAC, and others) are being used to predict properties of possible new DESs, a significant portion of publications related to these methods focus on correlating experimental data from type III DESs. [START_REF] González De Castilla | Thermodynamic and Transport Properties Modeling of Deep Eutectic Solvents: A Review on gE-Models, Equations of State, and Molecular Dynamics[END_REF] The theoretical studies of the structure and dynamics of DESs rely on all-atom (AA) molecular dynamics simulations. These simulations allow for great accuracy but are resource intensive, making proper sampling in larger or dynamically slower systems more cumbersome. As an alternative to AA models, coarse-grained (CG) models have found widespread applications. [START_REF] Merchant | Annual Reports in Computational Chemistry[END_REF][START_REF] Joshi | A review of advancements in coarse-grained molecular dynamics simulations[END_REF] By replacing groups of atoms into effective interaction sites, so-called CG beads, a large computational speedup is obtained at the cost of a somewhat reduced accuracy. [START_REF] Ingólfsson | The power of coarse graining in biomolecular simulations[END_REF] One of the most popular CG models is the Martini force field, [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF] in which 2-4 heavy atoms and associated hydrogens are united into CG beads.

Martini has been successfully employed in studies of various biologically relevant systems [START_REF] Melo | High-Throughput Simulations Reveal Membrane-Mediated Effects of Alcohols on MscL Gating[END_REF][START_REF] Moiset | Disaccharides Impact the Lateral Organization of Lipid Membranes[END_REF][START_REF] Ingólfsson | Lipid Organization of the Plasma Membrane[END_REF][START_REF] Souza | Protein-ligand binding with the coarsegrained Martini model[END_REF][START_REF] Dadsena | Ceramides bind VDAC2 to trigger mitochondrial apoptosis[END_REF] and more recently in material sciences. [START_REF] Vazquez-Salazar | Martini coarse-grained models of imidazolium-based ionic liquids: from nanostructural organization to liquid-liquid extraction[END_REF][START_REF] Alessandri | Resolving Donor-Acceptor Interfaces and Charge Carrier Energy Levels of Organic Semiconductors with Polar Side Chains[END_REF][START_REF] Piskorz | Nucleation Mechanisms of Self-Assembled Physisorbed Monolayers on Graphite[END_REF][START_REF] Alessandri | Bulk Heterojunction Morphologies with Atomistic Resolution from Coarse-Grain Solvent Evaporation Simulations[END_REF] One key reason for the high applicability of Martini has to do with the top-down process by which it was parameterized; [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF][START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] placing focus on matching thermodynamic parameters provides a level of transferability, en-suring at least some level of compatibility between newly parameterized models and existing models, as has been demonstrated in recent work with ILs. [START_REF] Vazquez-Salazar | Martini coarse-grained models of imidazolium-based ionic liquids: from nanostructural organization to liquid-liquid extraction[END_REF] This fact encouraged us to create the first Martini DES models, namely choline chloride -urea, tetrabutylammonium chloride -decanoic acid, and DL-menthol -acetic acid, based on the latest version of the force field (Martini 3.0 47 ). With the parameterization of these models, not only could we study the properties of these solvents themselves, but also their interactions with the wide array of pre-existing Martini models. The latter possibility opens, among other things, a completely new realm of temporal and spatial resolutions for studies of DES liquid-liquid extraction (LLE) processes with molecular dynamics, as will be demonstrated later in this document.

The paper is structured as follows: We first describe the parameterization and reasoning behind these models. This is followed with a comparison of the properties of these DESs against experimental data and AA simulations. Finally, we demonstrate the performance of the new models with test cases in which we evaluate the performance of our models against experimental data in more complex environments, namely the extraction of aromatic compounds from oil and water, and the self-assembly of surfactants.

Methods

Computational Details

All simulations were performed with GROMACS 2018.1 [START_REF] Abraham | GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF][START_REF] Páll | Tackling Exascale Software Challenges in Molecular Dynamics Simulations with GROMACS[END_REF][START_REF] Pronk | GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit[END_REF][START_REF] Hess | GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation[END_REF][START_REF] Van Der Spoel | GROMACS: Fast, flexible, and free[END_REF][START_REF] Lindahl | GROMACS 3.0: a package for molecular simulation and trajectory analysis[END_REF][START_REF] Berendsen | GROMACS: A message-passing parallel molecular dynamics implementation[END_REF] using the Bussi-Donadio-Parrinello thermostat (V-rescale). [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] Simulation pressures were controlled with the Berendsen barostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] during equilibration, and with the Parrinello-Rahman barostat [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] during production simulations. Constraints, where applied, were solved using LINCS, [START_REF] Hess | LINCS: A Parallel Linear Constraint Solver for Molecular Simulation[END_REF] with a LINCS order of 4.

All-atom simulations were run with a 2 fs timestep, while the CG simulations used a 20 fs timestep. Electrostatic interactions in AA and CG systems were computed using the smooth particle mesh Ewald (PME) method, [START_REF] Essmann | A smooth particle mesh Ewald method[END_REF] although all CG systems were also tested with the Reaction-Field method. All atom models were based on the CHARMM36 force field, [START_REF] Huang | CHARMM36m: an improved force field for folded and intrinsically disordered proteins[END_REF][START_REF] Best | Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone , and Side-Chain 1 and 2 Dihedral Angles[END_REF] released in March 2019 and ported for GROMACS by Bjelkmar et al. [START_REF] Bjelkmar | Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models[END_REF]

Simulation setup

AA reference simulations of individual DES components were performed by placing the solute in question inside a 3.5 nm 3 cubic box and solvating the system with TIP3P water. In case the solute had a net charge, a counterion was used to keep the system neutral. After energy minimization, using steepest descend algorithm in GROMACS, the system was equilibrated first in a NVT ensemble at ambient temperature for 1 ns, then in a NPT ensemble and ambient temperature and pressure for 1 ns, followed by a production simulation in the same ensemble for 300 ns. The routine for CG models consisted of the same energy minimization step, followed by four equilibration steps with timesteps of 1, 2, 10, and 20 fs respectively.

The equilibration simulations totalled to 80 ns, and were followed by a 1 µs production simulation. Both the equilibration and the production simulations were performed with a NPT ensemble. CG simulations of pure DESs were set up in two ways; by either randomly placing each molecule within the simulation volume as demonstrated in Figure 1, or by combining two separate phases formed in two separate simulations. Energy minimization and equilibration for every CG system followed the same routine described above. All LLE simulations were set up by randomly placing molecules in the simulation volume. The detailed compositions of the systems are given in Supplementary Information section 5.

A special case was made for testing the properties of solid tetrabutylammonium chloride (TBAC). A unit cell of the crystal structure was obtained from the Cambridge Crystallographic Data Centre (CSD Entry: XUBGIN) and packed in to a 5×5×5 crystal using Mercury. [START_REF] Macrae | Mercury 4.0: from visualization to analysis, design and prediction[END_REF][START_REF] Macrae | Mercury CSD 2.0 -new features for the visualization and investigation of crystal structures[END_REF][START_REF] Macrae | Mercury: visualization and analysis of crystal structures[END_REF][START_REF] Bruno | New software for searching the Cambridge Structural Database and visualizing crystal structures[END_REF][START_REF] Taylor | Rules governing the crystal packing of mono-and dialcohols[END_REF] The atomistic crystal structure was subsequently forward-mapped using Py-CGTOOL, [START_REF] Graham | PyCGTOOL: Automated Generation of Coarse-Grained Molecular Dynamics Models from Atomistic Trajectories[END_REF] and simulated with anisotropic pressure coupling both, as a pure crystal and as a crystal immersed in disordered TBAC. The system setup is illustrated in Figure 2.

Parameterization

Initial CG models presented in this work were parameterized using a more recent development version of Martini 3 open-beta. [START_REF] Souza | 3 -Open Beta-Release[END_REF] All systems presented here were tested again with the published parameters of Martini 3, [START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] and all reported values are from simulations with these parameters. Here we only give a cursory account of the parameterization process. A more detailed description of parameterizing small molecules for the Martini 3 force field can be found in the recent work of Alessandri et al. [START_REF] Alessandri | Martini 3 Coarse-Grained Force Field: Small Molecules[END_REF] CHARMM36 and GROMACS compatible AA starting structures and topologies were obtained from SwissParam tool, 71 and used as a basis for determining bonded interactions for the corresponding CG model. Being aware of the fact that SwissParam output is primarily meant for molecular docking and other applications in which high-level details are not the highest priority, 71 we compared the bonded parameters produced by SwissParam to those produced by Q-Force, [START_REF] Sami | Organic Materials in Silico: From force field development to predicting dielectric properties[END_REF] a QM based toolkit which derives bonded parameters through Hessian fitting. Comparison of the two predicted topologies was carried by simulating both models for 250 ns at the AA level, after which they were mapped to the corresponding CG structure using the centre-of-geometry (COG) mapping method. Comparing how the difference in AA bonded terms carries out to the final CG model, we found the average deviation between the bond lengths in the two CG models to be around 2.22%, while the angles deviated by 2.73%. Results of these comparisons are presented more thoroughly in Supplementary Information section 2. Based on these findings, we concluded that the possible errors present in the bonded terms predicted by SwissParam are at a tolerable level.

Following the COG mapping and assignment of bonded interactions between CG beads, the chemical nature of each moiety (represented by a bead) was evaluated and a relevant bead type was assigned to it. Fit of the bonded parameters (bond lengths, reference angles etc.) was improved by iteratively changing the parameters until a reasonable fit between the mapped AA models and the CG models was achieved. Finally, the shape and occupied volume of the CG model was semi-quantitatively compared to its atomistic counterpart via computation of the solvent accessible surface area (SASA) and the resulting Connollysurface. We used two sets of approximate values for atomistic radii; the default values distributed with most GROMACS versions, [START_REF] Bondi | van der Waals Volumes and Radii[END_REF] and a more recent set from Rowland et al. [START_REF] Rowland | Intermolecular Nonbonded Contact Distances in Organic Crystal Structures: Comparison with Distances Expected from van der Waals Radii[END_REF] These results are reported in the Supplementary Information section 3. Partial charges of choline, tetrabutylammonium and urea were derived using the CHelpG scheme, [START_REF] Breneman | Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis[END_REF] computed with Gaussian16 software. 76

Free energy estimates and partition coefficients

Assessing the correctness of non-bonded parameters, modelled via the Lennard-Jones and the Coulomb interactions, for individual CG models was performed by comparing experimentally and theoretically derived water-octanol partition coefficients (log P) to values computed for our models by thermodynamic integration (TI). The TI routine is described in full detail in the Supplementary Information section 4.

Thermal expansion coefficient

We estimated α, the thermal expansion coefficient of each DES studied, and compared the results against values reported in literature wherever available. α is given by:

α = - ( ∂ lnρ ∂ T ) p ( 1
)
where ρ is the density of the evaluated system, T is temperature, and p is pressure. In order to evaluate α, each DES was simulated three times as either dry (χ

H 2 O = 0 ) or wet (χ H 2 O = 0.1 )
within a temperature range of 288 -338 K. The range was scanned in intervals of 5 K. In the case of DL-menthol -acetic acid DES the upper limit of the temperature range was extended from 338 K to 353 K due to available experimental data of densities at higher temperatures.

Radial distribution functions

Liquid structures of CG DES systems were studied by comparing them against their atomistic counterparts by evaluating the radial distribution functions of key component pairs:

hydrogen bond acceptor (HBA) -HBD, HBA -water, and HBA -counterion. Bin volume and density were normalized, and the evaluation was performed on the centre-of-mass of each component, in order to obtain comparable results between the two resolutions.

Extraction efficiency

In order to compare the results from our LLE simulations to experiment, we computed the extraction efficiencies for each system. We calculated the extraction efficiency with:

f 0 -f 1 f 0 × 100% (2)
where f is the number density of the solute which is to be extracted. f 0 refers to the initial density of the solute within the solvent phase at the start of the process and f 1 to the density of the solute in the same phase after the process. We computed these values using the density profile of our system (typically along the z-axis), and by evaluating the number density of the target compound within the bulk of both phases, ignoring the density within the interface region.

Asphericity and shape estimations

The last test case presented in this study focused on sodium dodecyl sulfate micelles (SDS),

and their conformations when solvated in water or choline chloride -urea DES. The lengths of the micelles were computed from a distance matrix consisting of positions of all SDS molecules at every point in time. The radius of gyration and moments of inertia were computed using built-in GROMACS tools. Asphericity was computed using the scheme published by Dima et al., [START_REF] Dima | Asymmetry in the Shapes of Folded and Denatured States of Proteins[END_REF] which has been implemented in the MDAnalysis toolkit. [START_REF] Gowers | A Python Package for the Rapid Analysis of Molecular Dynamics[END_REF][START_REF] Michaud-Agrawal | A toolkit for the analysis of molecular dynamics simulations[END_REF]

Viscosity calculation

Shear viscosities of the simulated DESs were estimated with the Einstein relation, which is implemented in the gmx energy tool, available in GROMACS:

48 η = 1 2 V k B T lim t→∞ d dt ⟨( ∫ tt 0 +t t 0 P xz (t ′ )dt ′ ) 2 ⟩ t 0 (3)
In the equation above, V is the system volume and P xz is the average of the three off-diagonal elements of the pressure tensor.

Following the approach reported by Maginn et al., [START_REF] Maginn | Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0[END_REF] we evaluated the viscosity of each DES at eight equally spaced temperatures, ranging from 283, 293, ... , 353 K. A hundred replicas were run for each temperature, totaling 2400 separate simulations. Each system was equilibrated as described previously, and simulated in a NVT ensemble for 2.5 nanoseconds with a timestep of 5 fs.

Surface tension calculation

To assess the volatility of our DES models, and the correctness of our intermolecular interactions, we computed surface tension (γ) given by:

γ = L z 2 ⟨ P zz - P xx + P yy 2 ⟩ (4)
where L z is the length of the z-axis, P xx ,P yy , and P zz are the diagonal components of the pressure tensor, corresponding to the X, Y, and Z-axis, respectively. The factor 1/2 comes from the number of interfaces present in the simulation, which was two in our systems. Evaluation of the surface tension of each DES was done by equilibrating a system as described previously, and then elongating Z-axis such that over 66 % of the simulation volume consists of vacuum. For each DES, 20 replicas where simulated in 298 K for 200 nanoseconds, after which the surface tension was computed.

Results

Construction of DES topologies

Our study begun with the parameterization of five new CG models of hydrogen bond acceptors (HBAs) and donors (HBDs), allowing for the simulation of three DESs; the hydrophobic tetrabutylammonium chloride -decanoic acid DES, the hydrophilic choline chloride -urea DES, and the DL-menthol -acetic acid NADES. Acetic acid had been previously parameterized for Martini 3 by Souza et al. [START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] In addition, we parameterized molecules for the application of these DESs, namely benzothiophene, 2-methylthiophene, imidacloprid and dodecyl sulfate. CG models of all the components with their respective chemical structures are given in Figure 3. Full details of all newly parameterized models can be found in the Supplementary Information section 1.

Most of the topologies presented here can be derived simply by following the established rules, which can be found in the recent Martini 3 paper, [START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] or online at cgmartini.nl in the 'Tutorials' section. Since explicit hydrogen bonding is not possible at the CG resolution, these models require the use of the hydrogen bond donor and acceptor labels (d and a, respectively), which capture the hydrogen bonding propensities of the underlying chemical fragments in an effective manner. In the case of tetrabutylammonium cation, the standard mapping scheme resulted in the central nitrogen and first carbons of each four butyl chains being mapped in to a single regular Q-bead, carrying a positive unit charge. The remaining carbons on each butyl chain were mapped to SC3 beads, which in turn shielded the central Q-bead and severely hindered electrostatic interactions between the charged Q-bead and the counterion. To remedy this, we constructed four virtual sites in tetragonal arrangement around the central bead and divided the charge between these four sites, and the central bead. This arrangement was found to significantly improve the interaction between TBA + and Cl -, and allowed the central bead to be mapped to a more characteristic N2q-bead, due to the decreased partial charge on that bead. Further details are given in the Supplementary Following the general rules for mapping, urea was originally mapped as a single SP2 bead.

However, this would have made it almost indistinguishable from methanol, and completely indistinguishable from acetic acid. This was circumnavigated by placing two non-interacting U-beads were at 0.1 nm from the SP2 bead which were assigned opposite and equal partial charges in order to at least partially capture the highly polar nature of the compound.

Further details are given in the Supplementary Information section 1.5.

Validation of the DES models Density, thermal expansion coefficients, viscosity and surface tension

Densities of pure DESs were evaluated within a temperature range of 288, 293, 298... 338 K and compared against experimentally available data. Both, a completely dry, and water saturated systems were simulated in three replicas, separately for each DES. The resulting density curves as a function of temperature are given in section 6 of the Supplementary Information. Generally the Martini 3 models tend to overestimate the density by 5 to 10 %. The temperature dependency of the densities are captured reasonably well, as is evident from the thermal expansion coefficients (α) presented in Table 1. Predicted shear viscosities were systematically lower than those reported from experiment, [START_REF] Yadav | Densities and Viscosities of (Choline Chloride + Urea) Deep Eutectic Solvent and Its Aqueous Mixtures in Temperature Range 293[END_REF][START_REF] Van Osch | Hydrophobic deep eutectic solvents as water-immiscible extractants[END_REF][START_REF] Martins | Tunable Hydrophobic Eutectic Solvents Based on Terpenes and Monocarboxylic Acids[END_REF] but some general trends were captured. Further information is given in the Supplementary Information section 7. Predicted surface tensions match experimental values [START_REF] D'agostino | Molecular motion and ion diffusion in choline chloride based deep eutectic solvents studied by 1H pulsed field gradient NMR spectroscopy[END_REF][START_REF] Omar | Novel Nonanol-Based deep eutectic solvents: Thermophysical properties and their applications in Liquid-Liquid extraction and amino acid detection[END_REF][START_REF] Al-Akayleh | Therapeutic deep eutectic system of capric acid and menthol: Characterization and pharmaceutical application[END_REF] relativelty well, with all values being within 20 mN × m -1 from their experimental counterparts. Further information is given in the Supplementary Information section 8.

Partition coefficients

A key hallmark of the Martini model is reproduction of partition coefficients. To validate our DES models in this regard, water-octanol partition coefficients (log P) were computed using TI, and matched against experimental and predicted values with the exception of dodecyl sulfate for which no experimental log P values could be found. Results are shown graphically in Figure 4, and the corresponding data is given in Table 2. On average, the log P values of the Martini 3 models differed from experimental values by 0.2 log P units, corresponding roughly to half k B T . Choline-cation was the only significant outlier, differing by roughly one log P unit from the experimental value.

Liquid structure

Another key factor in assessing the correctness of our model was to compare the liquid structure within DESs against structures predicted by AA simulations. We placed our focus on radial distribution functions (RDFs) of the quaternary ammonium cations with respect to their counterions, water and the HBD. The results, as shown in Figure 5, indicate that the RDFs are in a good agreement with each other.

Figure 5: Radial distribution functions of TBA + -Cl -and Choline + -Cl -on both, coarsegrained and atomistic level.

Properties of solid TBAC

Most DES components remain in a solid state in ambient conditions, with some having melting points well above room temperature. Although the Martini model is primarily parameterized for modelling systems in the liquid state, we investigated whether or not a solid state can actually be reproduced with our DES topologies. As a test case, we selected a common DES component, TBAC, due to the availability of information about its crystalline structure and its common use in DES studies. Melting point of crystalline TBAC was estimated by simulating a single TBAC crystal immersed in an equal amount of fluid TBAC, in temperatures ranging from 250 to 315 K, with temperature increments of 5 K. From these simulations we expect to see the crystal phase growing below the melting temperature, and a fluid phase above the melting point. Each temperature was simulated three times for 300 ns, with the exception of 315 K, which results in a rather unstable system, due to the anisotropic pressure coupling. For each simulation, we computed the diffusion coefficient of both ions (TBA + and Cl -), since partial and complete melting of the system is easily observed from these values. While computing the MSD, no distinction was made between the compounds originally residing in the crystal and in the fluid phase. Results from this analysis are given in Figure 6. Experimentally the melting point has been estimated to be 310±1 K, [START_REF] Coker | Fusion properties of some ionic quaternary ammonium compounds[END_REF] and our model predicts the melting temperature to be around 300 to 305 K, at which point we notice sharp increases in the magnitude of the diffusion coefficients, from 10 -5 to 10 0 cm 2 /s, as well as in the overall cell volume (not shown). This increase signals the formation of a fluid phase as can also be seen in the graphical snapshots in Figure 6. In order to test the predictive capability of our DES models, we parameterized Martini 3 models for two target compounds: Benzothiophene (BTP) and 2-methylthiophene (MTP). These were solvated in octane, which has been used as a model-fuel in previous experimental studies alongside heptane, which is almost indistinguishable from octane at the coarse-grained level. We then proceeded to extract both, MTP and BTP from the model fuel by using a TBAC -acetic acid DES in 1:2 ratio. The octane-DES mass ratio was approximately 1:1. The choice of using AcAc as a HBD, as well as the octane-DES mass ratio, were based on the experimental work of Li et al., described in their 2016 study on the same topic, with the exception of using tetrabutylammonium bromide instead of TBAC. 106 The details of the simulation setup are given in Supplementary Information section 5. As shown in the graphical snapshot of the system (Figure 7), the octane-DES system forms a stable biphasic system starting from a homogeneously mixed initial state.

Both Our results suggest that the latter of the proposed mechanisms is more pronounced, since the extraction in our models is primarily driven by the favourable SC3 -TCn (n = 2, 4, 5, 6) interactions, which were designed to at least partially capture the CH-π interactions. In all likelihood the mechanism also relies on hydrogen bonding, which in Martini 3 is accounted for in the Lennard-Jones potential parameters.

Imidacloprid extraction from water

Neonicotinoids are a class of insecticides, previously lauded for their high efficacy and low risks. 108 Recent discoveries have pointed towards a unwanted effects in non-target organisms, especially in pollinators, 109 leading to a ban in some EU-countries and in the US. Although the search for new and alternative insecticides is well on its way, 110 neonicotinoids are still in use in many countries, 111 some of which are reporting frequent detection of neonicotinoids in drinking water. 112 Modern water treatment facilities use multi-step processes to remove potential unwanted residues from wastewater, but adsorption processes are still primarily used for removal of micropollutants, such as neonicotinoids. 113 Typically these processes use active carbon as an adsorbent, but the successful use of ILs has also been reported in cases where more specificity is required. We set up simulations using six equally spaced mole fractions for SDS, ranging from χ SDS = 0.0025 to χ SDS = 0.015, randomly placing them in a 1:2-ratio choline chloride -urea DES. Each simulation lasted 3 µs. Reference simulations were performed, in which same amounts of SDS were solvated in water for an equal time. For each system, we computed the length of the micelle, the radius of gyration and moments of inertia, and asphericity. We also report the number of micelles present in the system in the last 1 µs of the simulation. All numerical data is reported in the Supplementary Information section 9. The results show that our models reproduce the phenomenon reported by Arnold et al.: The formation of cylinder-like SDS micelles in the DES, and spherical SDS micelles in water. The behaviour in water and DES is illustrated in Figure 9. As a final test, one of the elongated micelles was removed from the DES phase and solvated with water in order to study whether it would remain stable in the aqueous phase. As is evident from the inset in Figure 9, the micelle rapidly disintegrates and begins to form spherical micelles. This finding is supported by computational work done by Tang et al.,118 who report that the bulk of the sodium is residing within 6 Å of the micelle when solvated in water.

Surfactant self-assembly in choline chloride -urea DES

Taken together, our results point to a replacement of sodium by choline cations when the micelles are solvated in choline chloride -urea DES. Given the much higher concentration of choline cations inside the DES with respect to sodium ions in the aqueous phase, the interactions between the charged SDS headgroups are much better screened, changing the effective curvature of the molecules toward less positive values and driving the transition to a less curved assembly.

Conclusions

We propose the first CG models for type III DESs and NADESs. Our models are capable of reproducing density, liquid structure and thermal expansion coefficients of pure DESs with relatively good agreement with experiment. We demonstrated the transferability of our models with three varied test cases, all of which were in good agreement with experiment. Compared to the existing atomistic force fields, our approach offers some important advantages: Coarse-graining allows for larger system sizes and longer simulations, reducing the need to use HPC facilities, and enables the study of more complex systems and processes, such as LLE, on more reasonable timescales. Furthermore, since the DES models are parameterized for Martini 3, following the same parameterization strategy, they are compatible with other Martini 3 models, and can be used in combination with biomolecular or polymeric systems. The availability of more systematic experimental data, in particular on partitioning of compounds between DESs and other solvents would facilitate the extension of our current DES topologies using the approach described in the paper.

As the library of validated DES compounds grows, Martini based simulations can be

Figure 1 :

 1 Figure 1: Example system setup for biphasic DES systems. (A) Random configuration starting structure. (B) Final self-assembled biphasic structure. The mixture shown in these images are composed of tetrabutylammonium chloride (blue and yellow, respectively), acetic acid (green), and octane (purple).

Figure 2 :

 2 Figure 2: System setup for TBAC crystal simulations. (A) Atomistic structure. (B) CG structure obtained by forward mapping using PyCGTOOL. (C) Crystalline TBAC (dark blue) immersed in disordered TBAC (cyan).

Figure 3 :

 3 Figure 3: DES compounds used in this study. (A) DES forming compounds, with the hydrogen bond acceptors (HBAs) on the top and hydrogen bond donors (HBDs) on the bottom. Compounds connected with a solid black line indicate DESs created in this study. A dashed line is drawn between TBA and AcAc, indicating their use in extraction of MTP and BTP. (B) Target compounds used in applications, which are described later in this paper. Red colour indicates a P-type bead, cyan corresponds to a C-type bead (with one exception, see SI section 1.8 for more information), blue is used for N-type beads, and yellow indicates a Q-type bead. Transparent beads are virtual sites or non-interacting U beads, which are described in more detail in the SI section 1.

Figure 4 :

 4 Figure 4: Computed, experimental, and predicted log P values of the 9 compounds modelled for this study.

Figure 6 :

 6 Figure 6: Melting point of crystalline TBAC. Diffusion coefficients of TBA + and Cl -are in blue and green, respectively. Partial melting begins at approximately 300 K, and the system reaches a fully molten state at 310 K

Figure 7 :

 7 Figure 7: Extraction of MTP and BTP from octane by TBAC -AcAc DES. (A) Evolution of the DES-octane mixture over time with snapshots of the system at various times, indicating the formation of a biphasic system. (B) Partial density graph of each component, along with a zoom in to the equilibrated system. In the partial density graph MTP, line drawn in green, can be seen favouring the DES phase, but still partially residing on the interface as well as in the octane phase. BTP, drawn in yellow, is almost completely extracted. (C) Components of the simulation, with their respective colours.

  114 Development of hydrophobic DESs has sparked interest in usage of DESs in wastewater treatment, with the hopes of finding a more recyclable and cheaper alternative for ILs.In 2017 Florindo et al.115 published their work on extracting acetamiprid, nitenpyram, thiamethoxam and imidacloprid from water using DL-menthol or TBAC based DESs. They demonstrate three successful extractions of all neonicotinoids from water with DL-menthol -carboxylic acid DESs, and report all TBAC based DESs leaching in to the water phase, contaminating it. We selected TBAC -DecA and DL-menthol -DecA as test cases, with the former reportedly failing due to TBAC leaching, and the latter being successful. For this particular test case we selected only imidacloprid (IMID) as a target. For DL-menthol -DecA systems, we set up four 2 µs simulations. Each simulation had χ IM ID =0.01 and 1:1 ratio between the DES components. The same was done for the TBAC -DecA systems, with the exception of 1:2 ratio between TBAC and DecA.Our results show similar behaviour to what Florindo et al. have reported: TBAC-DecADES extracts IMID, but simultaneously TBAC leaches in to the water phase. This is evident from the graphical snapshot as well as the partial density profiles, shown in Figure8.Menthol-DecA DES extracts IMID with an average efficiency of 75.2 ± 13 %, which is some-is partly due to the relatively high density of IMID within the interface region between the DES and water phase. The effect of the interface is much more pronounced in simulations than in experiment, due to the limited volume of our systems.

Figure 8 :

 8 Figure 8: Extraction of IMID from water by menthol -DecA DES. (A) Menthol -DecA DES extracting IMID from water. (B) Normalized partial density profiles of systems containing TBAC -DecA (left) or menthol -DecA (right). (C) TBAC -DecA DES in water, with TBA + leaching to the aqueous phase. (D) Compounds present in the simulation, with their corresponding colours.

Figure 9 :

 9 Figure 9: Asphericity of SDS micelles as a function of the mole fraction of SDS in DES.Corresponding values for water are given in the Supplementary Information section 9. Inset: An elongated micelle breaking to multiple spherical micelles when removed from the DES and solvated with water.

Table 1 :

 1 Thermal Expansion Coefficients of DESs studied in this work.

			Thermal Expansion Coefficients (α)
	DES	Water Martini 3 (×10 -3 K -1 ) Experiment (×10 -3 K -1 ) Reference
	TBAC -DecA	Dry	0.929 ± 0.004	1.040 ± 0.008	Ref. 81
	TBAC -DecA	Wet	0.931 ± 0.004	1.060 ± 0.005	Ref. 81
				0.696 ± 0.002	Ref. 82
	Menthol -AcAc Dry	1.671 ± 0.024	0.925 ± 0.010	Ref. 83
	Menthol -AcAc Wet	1.637 ± 0.026	0.886 ± 0.009	Ref. 83
	Chol. Cl -urea	Dry	0.617 ± 0.004	0.458 ± 0.006	Ref. 84
				0.533 ± 0.016	Ref. 85
				0.436 ± 0.010	Ref. 86
	Chol. Cl -urea	Wet	0.626 ± 0.002	0.531 ± 0.016	Ref. 85
				0.449 ± 0.001	Ref. 86

Table 2 :

 2 Log P values obtained with Martini 3, experimental values with their respective references, and predicted values with the software they were predicted with.

		Water -octanol partition coefficients
	Compound	Martini 3	Experimental Reference Predicted	Software
	TBA +	4.47 ± 0.06	4.51 †	Ref. 92	3.54	AlogPs
	DecA	4.52 ± 0.02	4.09	Ref. 93	4.10	XLogP3
	Menthol	2.73 ± 0.17	3.03	Ref. 94	3.40	XLogP3
	Choline +	-2.95 ± 0.04	-3.77	Ref. 95	-3.60	AlogPs
	Urea	-1.23 ± 0.04	-1.54	Ref. 96	-1.40	ChemAxon
	Benzothiophene	2.73 ± 0.03	3.12	Ref. 93	3.10	XLogP3
	2-methylthiophene 2.61 ± 0.05	2.33	Ref. 97	2.30	XLogP3
	Imidacloprid	1.07 ± 0.07	0.57	Ref. 98	0.41	KowWin
	Dodecyl sulfate -	4.63 ± 0.11	N/A		4.42	ChemAxon

† This value is obtained through QM computation, and thus it is not a true experimental value. However, the approach is built on higher level of theory than is typically used for predicting log P values. See reference nr. 92, p. 123 for more details. Simulation results are at 310 K; the temperature of the experimental measurements are not explicitly defined, and are assumed to be ambient.

  MTP and BTP are preferentially solvated in the DES phase, with BTP almost completely extracted, which is in line with experimental findings. BTP) individually, but at the same time. The exact method by which MTP and BTP are extracted is not fully understood. Li et al. propose hydrogen bonding as the primary mechanism for extraction of both BTP and thiophene, 100 whereas Wang et al. point towards CH-π interactions. 105

106

Based on 6 independent simulations, each lasting 2 µs, we obtained a 63.2 % ± 5.7 % extraction efficiency for MTP, and a 98.9 % ± 0.1 % extraction efficiency for BTP. Results and a comparison to experimentally derived extraction efficiencies are given in Table

3

.

We find our results in good agreement with the values reported in literature, in particular regarding the larger extraction coefficient of BTP with respect to MTP. There are a few considerations for reviewing this data: As this was just a 'proof of concept' type of test case, we did not screen various TBAC:AcAc ratios, nor DES:Oil ratios, which almost certainly would affect the extraction efficiency. In addition, we did not extract each component (MTP,

Table 3 :

 3 Extraction Efficiencies.

	HBA	HBD	Solute	Fuel	Extraction efficiency (%) Reference
	TBAC	Acetic Acid	BTP	Octane	98.9± 0.1	This work
	TBA-Br	Acetic Acid	BTP	Octane	81.75	Ref. 106
	TBAC Propanoic Acid	BTP	Octane	71 -99.47	Ref. 100
	TBAC	Acetic Acid	MTP	Octane	63.2± 5.7	This work
	TBA-Br	Acetic Acid	Thiophene Octane	55.8	Ref. 106
	TBAC	Malonic Acid	MTP	Octane	50.05	Ref. 107

  (5) Zhekenov, T.; Toksanbayev, N.; Kazakbayeva, Z.; Shah, D.; Mjalli, F. S. Formation of type III Deep Eutectic Solvents and effect of water on their intermolecular interactions. Fluid Phase Equilibria 2017, 441, 43-48, 2016 Asian Thermophysical Properties Conference. (6) Paiva, A.; Matias, A. A.; Duarte, A. R. C. How do we drive deep eutectic systems towards an industrial reality? Current Opinion in Green and Sustainable Chemistry 2018, 11, 81-85, Pharmaceuticals / Green Processes and Technologies. (7) Gill, I.; Vulfson, E. Enzymic catalysis in heterogeneous eutectic mixtures of substrates. Journal of the American Chemical Society 2004, 126, 9142-9147, PMID: 15264850. (11) Dai, Y.; Witkamp, G.-J.; Verpoorte, R.; Choi, Y. H. Tailoring properties of natural deep eutectic solvents with water to facilitate their applications. Food Chemistry 2015, 187, 14-19. (12) Abo-Hamad, A.; Hayyan, M.; Saadi, M. A. H. A.; Hashim, M. A. Potential applications of deep eutectic solvents in nanotechnology. Chemical Engineering Journal 2015, 273, Zhang, Q.; De Oliveira Vigier, K.; Royer, S.; Jérôme, F. Deep eutectic solvents: syntheses, properties and applications. Chem. Soc. Rev. 2012, 41, 7108-7146. (14) Carriazo, D.; Serrano, M. C.; Gutiérrez, M. C.; Ferrer, M. L.; del Monte, F. Deepeutectic solvents playing multiple roles in the synthesis of polymers and related materials. Chem. Soc. Rev. 2012, 41, 4996-5014. (15) Perna, F. M.; Vitale, P.; Capriati, V. Deep eutectic solvents and their applications as green solvents. Current Opinion in Green and Sustainable Chemistry 2020, 21, 27-33, New Synthetic methods Selected papers from the 4th Green and Sustainable Chemistry Conference. (16) Xu, G.-C.; Ding, J.-C.; Han, R.-Z.; Dong, J.-J.; Ni, Y. Enhancing cellulose accessibility of corn stover by deep eutectic solvent pretreatment for butanol fermentation. Recyclable deep eutectic solvent for the production of cationic nanocelluloses. Carbohydrate Polymers 2018, 199, 219-227.(19) Haghbakhsh, R.; Raeissi, S.; Duarte, A. R. C. Group contribution and atomic contribution models for the prediction of various physical properties of deep eutectic solvents. Scientific Reports 2021, 11, 6684. (20) Alkhatib, I. I.; Bahamon, D.; Llovell, F.; Abu-Zahra, M. R.; Vega, L. F. Perspectives and guidelines on thermodynamic modelling of deep eutectic solvents. Journal of Molecular Liquids 2020, 298, 112183. For Table of Contents Use Only. Synopsis: New CG models for DESs allow for larger scale simulations which offer more relevance for the development and testing of new sustainable solvents.

	(13) Bioresource Technology 2016, 203, 364-369.
	551-567.
	29

Trends in Biotechnology 1994, 12, 118-122. (8) Lopez-Fandino, R.; Gill, I.; Vulfson, E. N. Protease-catalyzed synthesis of oligopeptides in heterogenous substrate mixtures. Biotechnology and Bioengineering 1994, 43, 1024-1030. (9) Abbott, A. P.; Capper, G.; Davies, D. L.; Rasheed, R. K.; Tambyrajah, V. Novel solvent properties of choline chloride/urea mixtures. Chem. Commun. 2003, 70-71. (10) Abbott, A. P.; Boothby, D.; Capper, G.; Davies, D. L.; Rasheed, R. K. Deep Eutectic Solvents Formed between Choline Chloride and Carboxylic Acids: Versatile Alternatives to Ionic Liquids. (17) Lapena, D.; Errazquin, D.; Lomba, L.; Lafuente, C.; Beatriz, G. Ecotoxicity and biodegradability of pure and aqueous mixtures of deep eutectic solvents: glyceline, ethaline, and reline. Environmental Science and Pollution Research 2021, 28, 8812-8821. (18) Li, P.; Sirviö, J. A.; Asante, B.; Liimatainen, H.

Acknowledgement

PAV would like to thank R. Alessandri, F. Grünewald and A.H. de Vries for many helpful discussions. The MD simulations were carried out on the Dutch national e-infrastructure with the support of SURFsara.

Supporting Information Available

Model data and parameters, comparison of bonded terms, SASA numerical data, TI routine, system compositions, densities, SDS micelle numerical data, viscosity data.


Orlando Villegas

Germain Salvato Vallverdu

Brice Bouyssière

Socrates Acevedo
email: socrates.acevedo@gmail.com

Jimmy Castillo

Brice Bouyssiere

Sócrates Acevedo

Isabelle Baraille

Cancellation of dipole moment of models of asphaltene aggregates as a mean for their dispersion in toluene and THF calculated using molecular dynamics

Cancellation of dipole moment of models of asphaltene aggregates as a mean for their dispersion in toluene and

THF calculated using molecular dynamics

Introduction

In view of the deleterious effects that asphaltene aggregation causes in the different stages of crude oil processing, a very large number of research reports dealing with this subject could be found in the literature and a very small number of these could be found in the Reference part below. [START_REF] Headen | Simulation of Asphaltene Aggregation through Molecular Dynamics: Insights and Limitations[END_REF][START_REF] Silva | 11 Asphaltene aggregation studied by molecular dynamics simulations: role of the molecular architecture and solvents on the supramolecular or colloidal behavior[END_REF][START_REF] Jiang | Molecular Mechanisms of Suppressing Asphaltene Aggregation and Flocculation by Dodecylbenzenesulfonic Acid Probed by Molecular Dynamics Simulations[END_REF][START_REF] Headen | Predicting Asphaltene Aggregate Structure from Molecular Dynamics Simulation: Comparison to Neutron Total Scattering Data[END_REF][START_REF] Headen | Evidence for Asphaltene Nanoaggregation in Toluene and Heptane from Molecular Dynamics Simulations[END_REF][START_REF] Villegas | Molecular Cartography of A1 and A2 Asphaltene Subfractions from Classical Molecular Dynamics Simulations[END_REF][START_REF] Silva | Molecular Dynamics Study of Nanoaggregation in Asphaltene Mixtures: Effects of the N, O, and S Heteroatoms[END_REF][START_REF] Chmiela | Towards exact molecular dynamics simulations with machinie-learned force fields[END_REF][START_REF] Schuler | Descriptive accounts of thermodynamic and colloidal models of asphaltene flocculation[END_REF] No doubts such interest has provided a large knowledge on the theme and controversies about basic concepts are common. For instance, whereas many conjectures could be proposed to account for aggregate formation of different sizes there is no consensus regarding the stability of these aggregates in solution. In this work we used a molecular dynamics method and polycyclic aromatic models containing heteroatom to simulate both aggregate formation and stability. We found that our simulation is coherent with formation of aggregates with dipole moments (µ) equal or close to zero. Moreover, we also found that µ -→ 0 when the number of molecules within the aggregate tends to a large number.

The study of asphaltenes has always been involved in some difficulties due to their complex constitution. Because of this, his studies on solutions are debated in two main theories, the colloidal solution theory or colloidal model, 10 this theory helps to explain the presence of small, and dense nanoaggregates in solution. It has also been the pillar of many structural models. [11][12][13] The phase theory for multicomponent mixtures or thermodynamic model is also involved in the behavior of asphaltenes in solution, 14,15 especially in studies concerning the solubility parameter. [16][17][18] This shows the dual nature of asphaltenes depending on the concentration present, partially dissolved, partially colloidal. [19][START_REF] Balestrin | Recent Developments on the Elucidation of Colloidal Aspects of Asphaltenes and Their Relevance to Oilfield Problems[END_REF][21] Molecular dynamics simulations provide us with an ideal technique to study the behavior of a system at nanometer scales, where intermolecular interactions are present. These simulations allow us the sampling of a system in time by generating a set of configurations where the conditions of temperature, pressure and species are constant. This is helpful in observing the behavior of asphaltene aggregates in solution. 10,14 However, due to the large family of compounds present, it is difficult to devise structural models. Recent work, using Atomic Force Microscopy (AFM) has isolate real asphaltene molecules 9 which resulted to be polycyclic aromatic compounds (PAC) largely similar to empirical aromatic models used to represent the behavior of asphaltenes in solution. [START_REF] Villegas | Molecular Cartography of A1 and A2 Asphaltene Subfractions from Classical Molecular Dynamics Simulations[END_REF] Using the below described molecular dynamics methods and appropriate molecular models we calculated here that cancellation of dipole moments of aggregates is a principal factor in the stability of asphaltene aggregates containing many molecules.

Methods

Two different models, type A1 and type A2 were used for asphaltene simulations. Figure 1 shows the 3D molecular models employed for A1, A2 and aggregate composed by 50% of each one. The molecular dynamics calculations were performed in toluene and THF at room conditions. 80 asphaltenes molecular models were randomly introduced in an 8x8x8 nm3 toluene or THF box which was equilibrated during 60 ns, before performing the molecular dynamics. Once the molecular dynamics of the system were completed the aggregates formed in the system were analyzed and characterized. All simulations were carried out using the GROMOS96 force fields with 53A6 parameter set, 11 using the GROMACS 2019.3 computer package [START_REF] Headen | Evidence for Asphaltene Nanoaggregation in Toluene and Heptane from Molecular Dynamics Simulations[END_REF]11 . The main feature of this force field is its united-atom (UA) parameter which allows an increase in the size of system without losing computer time. Other characteristics Figure 1: Structural representation of models A12 and A22 and aggregate with a 50 % mixture of both.The yellow spheres are S, the blue N and the red OH groups Figure 1 shown the molecular structure base of the asphaltenes molecules used in this study for the subfractions A1 and A2. All models have the same basic structure, consisting of an aromatic core, with a general island-like architecture. These are based on those proposed by Acevedo, [START_REF] Acevedo | Structural Analysis of Soluble and Insoluble Fractions of Asphaltenes Isolated Using the PNP Method. Relation between Asphaltene Structure and Solubility[END_REF] shown in Figure 1. The main difference between the A1 and A2 structures is the exchange of a side chain (n-octyl) for two alkyl cycles, which produces a decrease in the H/C ratio and increases the DBE producing a more rigid structure for A1 and a more flexible one for A2. [START_REF] Acevedo | Trapping of Paraffin and Other Compounds by Asphaltenes Detected by Laser Desorption Ionization-Time of Flight Mass Spectrometry (LDI-TOF MS): Role of A1 and A2 Asphaltene Fractions in This Trapping[END_REF] The presence of a sulfur atom, nitrogen atom, hydroxyl group on the aromatic nucleus and the addition of a carboxyl group on the side chain were considered.

Charge of atoms for all atoms in each aggregate were calculated using the semiempirical AM1 program. [START_REF] Jakalian | Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation[END_REF][START_REF] Thompson | Parameterization of charge model 3 for AM1, PM3, BLYP, and B3LYP[END_REF][START_REF] Jakalian | Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: I. Method[END_REF] These charges were then used for dipole moment calculation. These DM Figure 2 shows the dipole moments changes of the A12 aggregates relative to the number of molecules. Calculations were done for both toluene and tetrahydrofurane (THF). As depicted, when the number of model molecules in the aggregate increases their dipole moments decrease. In most cases the decrease is very fast. The dipole moments reach values close to zero for a number of molecules close to 45. In general, the dipole moment curve drops to a sort of plateau close to zero both for toluene and THF.

Discusion

Stability of aggregates in toluene, required for its successful dispersion has been calculated using molecular dynamics method, employing both toluene and THF. For this, cancellation of the dipole moment of aggregates is required. The dipole moment cancellation (µ = 0) is produced by randomly orientation of molecules within the aggregate. This cancellation would be the largest while the largest would be the number of calculated molecules inside the aggregate. In other terms: Dispersion stable µ = 0; when N is a large number N is the number of molecules in aggregate.

To our knowledge, in the literature there is no mention of the above condition (µ = 0) to account for dispersion stability of asphaltene. We found this to be the case in 22 experiments similar to the one described above for the A12 case. Although this is a welcome finding, more research is highly recommended before such a concept gain more consensuses. For instance what would be the results in the case of mixtures of molecular models.

According to the above results asphaltenes could be considered as self-assembling materials where the usual several physical factors such as hydrogen bonding, dispersion and polar effects promote aggregate formation and growing (increase N Agg ) and simultaneously such increase in promotes cancellation or annulment of dipole moment thus suppressing aggregate formation. This compromise between these opposite effects (N Agg ) would be responsible for the stability of the dispersion.

Conclusions

Calculation of dipole moments of model aggregates in solutions, both in toluene and THF, using the above described molecular dynamics technique constitutes a valid approach to understand both formation and stability of these solutions. In view of them in being a stable dispersion of large and small aggregates, the proposal regarding dipole moment cancellation, linked to intermolecular attraction of molecules within the aggregates completes the required steps to obtain the corresponding dispersion.

Figure 2 :

 2 Figure 2: Reduction of the dipole moment with the increase of the number of molecules in the aggregate according to the simulation. Here we observe similar results for models A12 and A22 in toluene and THF.

Figure 3

 3 Figure 3 shows various arrays of aggregates for the A12 model. Dipole moments of these arrays and corresponding number of molecules are listed in Table 1. As the number of molecules increases, the dipole moment decreases. Arrows indicate the direction of the dipole moment corresponding to molecular models. Random array of these leads to cancellation of the dipole moment of the aggregate.

Figure 3 :

 3 Figure 3: Array of molecular models corresponding to a single molecule and various aggregates for A12 model.

Table 1 :

 1 DM as a function of N DM as a function of N N It should be noted that dependence of dipole moment with N in THF is similar to the one observed in toluene see below Figure 4.

Figure 4 :

 4 Figure 4: Mixtures of the model shown in THF solvent.

Figure 4

 4 Figure 4 shown mixtures of the model in THF solvent. As was the case for single models shown in Figures 4, the dipole moment of aggregates also tends to zero when the number of molecules increases. In the Figure 5 is presented a sketch of the array of molecular models corresponding to a single molecule and various aggregates. As shown in Table 1, as the number of molecules increases, the dipole moment decreases. Arrows indicate the direction of the dipole moment corresponding to molecular models. Random array of these leads to cancellation of the dipole moment of the aggregate

Figure 5 :Figure 6 :

 56 Figure 5: Section of Molecular display corresponding to mixture A12 (in brown) +A22 (in pink. Solvent (not shown) is toluene

lar

  Structures of Asphaltenes by Atomic Force Microscopy. Journal of the American Chemical Society 2015, 137, 9870-9876. (10) Leontaritis, K. J.; Mansoori, G. A. Asphaltene Flocculation During Oil Production and Processing: A Thermodynamic Colloidal Model. SPE International Symposium on Oilfield Chemistry. 1987; p 16258. (11) Mullins, O. C. The Modified Yen Model. Energy & Fuels 2010, 24, 2179-2207. (12) Gray, M. R.; Tykwinski, R. R.; Stryker, J. M.; Tan, X. Supramolecular Assembly Model for Aggregation of Petroleum Asphaltenes. Energy & Fuels 2011, 25, 3125-3134. (13) Acevedo, S.; Castillo, J.; Vargas, V.; Castro, A.; Delgado, O.; Cortés, F. B.; Franco, C. A.; Bouyssiere, B. Suppression of Phase Separation as a Hypothesis to Account for Nuclei or Nanoaggregate Formation by Asphaltenes in Toluene. Energy & Fuels 2018, 32, 6669-6677. (14) Painter, P.; Veytsman, B.; Youtcheff, J. Asphaltene Aggregation and Solubility. Energy & Fuels 2015, 29, 2120-2133. (15) Andersen, S. I.; Speight, J. G. Thermodynamic models for asphaltene solubility and precipitation. Journal of Petroleum Science and Engineering 1999, 22, 53-66. (16) Acevedo, S.; Castro, A.; Vásquez, E.; Marcano, F.; Ranaudo, M. A. Investigation of Physical Chemistry Properties of Asphaltenes Using Solubility Parameters of Asphaltenes and Their Fractions A1 and A2. Energy & Fuels 2010, 24, 5921-5933. (17) Rogel, E. Theoretical Estimation of the Solubility Parameter Distributions of Asphaltenes, Resins, and Oils from Crude Oils and Related Materials. Energy & Fuels 1997, 11, 920-925. (18) Redelius, P. Bitumen Solubility Model Using Hansen Solubility Parameter. Energy & Fuels 2004, 18, 1087-1092.


Icaro V Soares
email: icaro.soares@ieee.org

Denys Nikolayev
email: denys.nikolayev@deniq.com

Stochastic Analysis of WPT Efficiency due to Location Uncertainty of mm-Sized Deep-Implanted Pacemakers

Keywords: bioelectronics, efficiency bounds, implantable devices, pacemaker, wireless power transfer

The lossy and heterogeneous nature of the body tissues and the significant wave-impedance contrast as the electromagnetic wave propagates inside the body limits the fundamentally achievable efficiency of far-and midfield wireless charging of deep-body implantable devices. In addition, fluctuations in the implant position cause a further variation in the maximum power transfer efficiency. In this work, the analyzed problem consists of a wireless charged deep-implanted pacemaker considering an electric and a magnetic radiation source and a uniform distribution of pacemaker's positions inside the heart. The results show that the optimal frequency remains the same despite the position fluctuations. However, the efficiency bounds are strongly related to the nature of the source. Therefore, the proposed model and the obtained results serve as a gauge for designing implantable wireless power transfer systems and predicting variations in their efficiency.

I. INTRODUCTION

Implantable bioelectronics emerges as an outstanding advance in medicine, enabling innovative approaches from diagnostics to treatments. Furthermore, the application of wireless power transfer techniques to these devices extends their potential applications and enables significant miniaturization of these implants. For instance, inductive wireless charging techniques have been used for powering sub-cutaneous pacemakers [START_REF] Liu | An effective sandwiched wireless power transfer system for charging implantable cardiac pacemaker[END_REF], [START_REF] Campi | Wireless power transfer charging system for aimds and pacemakers[END_REF], whereas powering the ones deep-implanted in the heart requires the application of far-and midfield approaches [START_REF] Ho | Wireless power transfer to deep-tissue microimplants[END_REF], [START_REF] Le-Huu | Bipolar spiral midfield wireless power transfer for cardiac implants application[END_REF].

However, further analysis of the electromagnetic wave energy transfer between the on-body transmitter and the implantable receiver is required to establish the efficiency limits and ensure compliance with the safety requirements. In order to determine the efficiency bounds on the power radiated by implantable devices, different approaches have been proposed in the literature. Most of them consist of stratified models with simplified geometry, in general, planar [START_REF] Poon | Optimal frequency for wireless power transmission into dispersive tissue[END_REF] or spherical [START_REF] Skrivervik | Fundamental limits for implanted antennas: maximum power density reaching free space[END_REF]- [START_REF] Nikolayev | Optimal radiation of body-implanted[END_REF]. These models predict the maximum theoretical efficiency This study was supported in part by ANR "MedWave" project and in part by the Région Bretagne through the project SAD "EM-NEURO". and the optimal frequency and quantitatively describe the loss mechanism in these systems. The simplified models lead to more generalized results; however, they cannot represent the complexities of the body structures. For this reason, some studies employ anatomical models [START_REF] Nikolayev | Electromagnetic radiation efficiency of body-implanted devices[END_REF]- [START_REF] Soares | Wireless powering efficiency assessment for deep-body implantable devices[END_REF] in their analysis when determining the maximum efficiency in a particular application is envisaged.

In this work, the problem of wirelessly charging a deepimplanted pacemaker is modeled based on a two-dimensional anatomical model of the human pectoral. Based on this model, the maximum achievable efficiency and the optimal frequency are calculated considering an electric and a magnetic transmitter, and their performance is compared. Besides, this model is used to evaluate the sensitivity of these efficiency bounds to the unceasing fluctuations in the implant positions caused by the heart movements. Therefore, the results described in this work can help design wireless powered deep-implanted pacemakers and predict their charging efficiency.

II. PROBLEM FORMULATION

In order to analyze the efficiency bounds on wireless powering of a pacemaker implanted directly in the human myocardium, the problem is geometrically represented through a two-dimensional cross-section of the pectoral region, as shown in Fig. 1(a). In this model, the closed regions representing the human organs are filled by dielectrics with equivalent complex permittivities as established in [START_REF] Gabriel | The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues[END_REF]. Furthermore, the transmitter and receiver are assumed to be perfectly aligned along the z-(out-of-plane) axis. Therefore, the translational invariance is justified, and the problem can be treated as twodimensional.

As indicated in Fig. 1(a), the receiver is represented as a lossless circular region with a radius R = 1 mm, filled by a dielectric with its real permittivity matched to the heart. In addition, the on-body transmitter is considered to surround the entire body surface representing the best-case scenario and therefore providing the fundamental efficiency bounds. In order to minimize the losses due to reflection at the interface between the body and external medium, the latter is filled by the same dielectric constant as the skin and zero conductivity.

Considering that the media composing the organs are linear and disregarding polarization mismatching, the reciprocal problem can be analyzed [START_REF] Jackson | Classical electrodynamics[END_REF]. Therefore, for calculation purposes, the radiation source is considered inside the circular domain Ω s , whereas the received power is evaluated over the on-body surface Ω p . In this way, based on the time-harmonic electric E and magnetic H fields, the transmitted P t , received P r , and dissipated power P d can be calculated through the integrals (1a)-(1c) below:

P t = Ωs 1 2 E × H * • ds, ( 1a
) P r = Ωp 1 2 E × H * • ds, (1b)
P d = 1 2 Ωp σ |E| 2 dv. ( 1c
)
So that the energy is conserved as follows [START_REF] Gabriel | The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues[END_REF]:

P t = P r + P + i2ω(W m -W e ), (2)
where W m and W e are the time-averaged energies stored in the magnetic and electric fields, respectively.

From the real parts of the transmitted ℜ(P t ) and received ℜ(P r ) powers, the wireless power transfer efficiency η can be defined as:

η ≡ ℜ(P r ) ℜ(P t ) . (3)
Disregarding the impedance mismatch between the antennas and the electronic circuits connected to them, as well as other system-level losses in these same circuits, the losses in an implantable wireless power transfer system are mainly due to the wave attenuation in the lossy body tissues. In addition, reflections at the interface body-transmitter, the heart-receiver, and between the inner organs reduce the power delivered to the receiver. While the reflection at the two first interfaces can be mitigated by adequately designing the transmitter and receiver (e.g., by employing high-permittivity and low-loss substrates [START_REF] Nikolayev | Dielectric-loaded conformal microstrip antennas for versatile in-body applications[END_REF] matched to the wave impedance at the skin and the heart, respectively), the losses due to attenuation and reflection at the transition between organs are unavoidable. Therefore, the way the problem is defined in this work assures that the calculated efficiency is limited only by the inexorable loss sources, so the following results characterize the efficiency bounds in the wireless charging efficiency for deep-implanted pacemakers.

III. METHOD OF ANALYSIS

The electromagnetic field and integrals (2a) and (2b) were numerically evaluated using the 2D full-wave simulation module implemented in the simulation suite COMSOL Multiphysics with the geometric representation shown in Fig. 1(a) and described in the previous section.

This work considers two radiation sources, an electric and a magnetic, as shown in Figures 1(b once in this direction, the radiation pattern is parallel to the chest boundary, in which the on-body region is the closest to the source. Consequently, this orientation maximizes the integral (1b).

Finally, in order to evaluate the effect of the implant positioning on the efficiency bounds, this simulation was performed with 200 uniformly distributed random positions inside the heart, such as represented in Fig. 1(a) by the gray circles.

IV. RESULTS

First, it is important to validate the assumptions that the reflections at the interface between the heart-receiver and the body-transmitter can be disregarded in this proposed model. The first case is justified as the receiver is located inside the heart, and its size is much smaller than the heart size. However, once the on-body transmitter occupies a larger surface, the reflections at the interface body and external medium are more significant. Besides, in this work, the external medium is filled by a lossless dielectric perfectly matched to the wave impedance in the skin. It can also be argued that the other tissues will change the effective permittivity, affecting the wave impedance inside the body.

For this validation, the reflection coefficient at the transition between two media is defined as:

Γ ≡ Z 1 -Z 2 Z 1 + Z 2 , (4)
where Z 1 and Z 2 are respectively the wave impedance in the medium 1 (inside the body) and 2 (external medium), which can be calculated as Z n = (jωµ 0 )/(σ n + jωε n ), being σ n the conductivity and ε n the electric permittivity of the medium n.

The absolute value for the reflection coefficient over a frequency span of 10 MHz to 10 GHz is presented in Fig. 2 for a plane wave propagating at the interface between two different media, considering the wave impedance in medium 1 the same as in the skin Z skin and also in the effective medium formed by the nine tissues Z effective . In contrast, the medium 2 was set as air Z air or the lossless and skin-like dielectric Z lossless skin .

As Fig. 2 demonstrates, there is a significant wave impedance contrast when the external medium is considered air, leading to a high reflection coefficient. Even though it decreases with the frequency, more than 70% of the transmitter power is reflected over the entire analyzed range. On the other hand, this reflection is considerably reduced when the external medium is matched to the body's wave impedance. For frequencies below 50 MHz, the electromagnetic properties of the other tissues affect the effective permittivity. Therefore, the reflection coefficient can be even twice as high if only the intrinsic impedance in this skin is considered. However, as the frequency increases, matching the wave impedance in the external medium to the skin leads to a numerically equivalent reflection coefficient matching with the effective medium. It is important to point out that the first approach is used once from the practical perspective is much easier to know the wave impedance only in the skin than experimentally characterize an effective medium. Finally, the results in Fig. 2 show that the assumptions proposed in this work are valid over 100 MHz and, in the range between 100 MHz and 5 GHz, the reflected power is lower than 20%. For this reason, this frequency span will be considered for the following analysis. After calculating the received power for each one of the 200 pacemaker positions indicated in Fig. 1(a), the efficiencies considering an electric and magnetic source are shown in Fig. 3(a) and 3(b), respectively. Overall, it can be noticed that a change in the receiver position inside the heart leads to a variation of the efficiency value inside a range, but the optimal frequency remains unaffected. Even though previous works have shown that the frequency in which the efficiency peaks reduces as the implantation depth increases [START_REF] Soares | Physical bounds on implant powering efficiency using body-conformal WPT systems[END_REF], the results presented in Fig. 3 are justified once for all positions, the receiver is confined in the same medium, and the implantation depth varies within a small range. Furthermore, the verified optimal frequency stability assures that these variations in the implant position cannot detune the system.

The nature of the source, in turn, substantially affects the efficiency bounds. First, for an electric source, the efficiency increases with the frequency until reaching its maximum at 1.097 GHz; after that, a sharp decrease is verified. At this optimum frequency, the average efficiency is equal to 1.85 × 10 -3 , but depending on the position, it goes from 9.49 × 10 -4 to 5.21×10 -3 . On the other hand, the magnetic source leads to maximum efficiency at low frequencies, and it decreases as the frequency increases. For comparison purposes, at 1.097 GHz, the average efficiency for a magnetic source is 1.03 × 10 -2 , ranging from 2.91 × 10 -2 to 5.05 × 10 -3 with the variation in the pacemaker's position.

The reflections at the interface between the body and the receiver or transmitter can be disregarded; also, the inner reflections equally affect the performance of both sources. Therefore, the factor that favors the most the magnetic source over the electric one is the fact that the body structures have low interaction with the magnetic field (note that the opposite trend can be observed for the case of implantable radiating structures [START_REF] Nikolayev | Optimal radiation of body-implanted[END_REF]). In contrast, the electromagnetic field is significantly attenuated in the tissues. Such attenuation increases with the frequency, but at lower frequencies, the near field losses deteriorate the η E . Consequently, a tradeoff between these two loss causes is reached at the optimal frequency [START_REF] Soares | Wireless powering efficiency assessment for deep-body implantable devices[END_REF].

V. CONCLUSION

In this work, theoretical bounds on the wireless charging efficiency of a deep-implanted pacemaker are evaluated as a function of the frequency and the implant position. First, the efficiency analysis over the frequency allows one to quantify the effect of losses on the efficiency and determine the optimal operating frequency in each case. In addition, by carrying out this analysis for many randomly distributed implant positions, it is possible to determine the variation in the maximum achievable efficiency due to fluctuations in the implantation depth.

Therefore, a model based on the two-dimensional crosssection of the human pectoral composed of nine different tissues is analyzed. Furthermore, this model is proposed so that all loss sources that could be mitigated in practical implementations are disregarded; thus, only the intrinsic losses are taken into account. In this way, the calculated efficiency corresponds to the maximum theoretical efficiency achieved in this given application.

The results obtained through this proposed model show that fluctuations in the implant position are not able to detune the system once the maximum efficiency is achieved at the same frequency for all positions. It is also shown that a magnetic radiation source monotonically decreases with the frequency, whereas an electric transmitter presents an optimum frequency at 1.097 GHz. This behavior is because there is practically no interaction between the magnetic field and the body tissues while strongly attenuating the electric field. The near field coupling reduces efficiency at lower frequencies for an electric source while the wave is almost completely attenuated at higher frequencies. Therefore, the efficiency reaches a peak at the mid-field range.

Even though this model theoretically leads to the maximum efficiency levels that could be attained, it also provides information about what could be done in a physical implementation. First, as the source ideally covers the entire body surface, it is reasonable to assume that conformal antenna array transmitters will lead to higher efficiency than single-antenna configurations [START_REF] Soares | Wireless powering efficiency assessment for deep-body implantable devices[END_REF]. Additionally, wavefront-shaping and focusing techniques will also increase the power delivered to the implanted receiver. For instance, the physical realization of a wirelessly powered pacemaker with a focused array transmitter proposed in [START_REF] Ho | Wireless power transfer to deep-tissue microimplants[END_REF] exhibited an efficiency of 3.90 × 10 -4 at the frequency of 1.6 GHz, which corroborates with the theoretical average levels presented in this paper.

Finally, it should be stated that a specific problem of a wirelessly charged pacemaker is investigated in this work; however, the proposed model can also be applied to analyze wireless power transfer efficiency in different implantable devices applications. Furthermore, the obtained results can be used as a guideline to design miniaturized deep-implanted pacemakers and predict wireless charging efficiency fluctuations.

Fig. 1 .

 1 Fig. 1. Proposed formulation for evaluating the efficiency bounds of wirelessly charging a deep-tissue pacemaker: (a) Anatomical two-dimensional model of a human pectoral region comprising nine lossy tissues. The transmitter is an on-body source that covers the entire body surface, whereas the receiver is a pacemaker implanted inside the myocardium. The gray circles represent the evaluated implanted positions. The sources are modeled as an (b) electric and (c) magnetic point dipole with moment vectors n E p and n H p , respectively. Their corresponding electric E (continuous blue line) and magnetic H (dashed pink line) fields are also indicated.

3 Fig. 3 .

 33 Fig. 3. Wireless power transfer efficiency for an (a) electric source and (b) magnetic source. The continuous blue line is the efficiency averaged (avg η) over the 200 evaluated receiver positions. The upper and lower limits of the blue area represent the efficiency at the position that leads to the maximum (max η) and minimum (min η) values. The red bar indicates the variation in the efficiency at the optimal frequency.

  Fig.2. Reflection coefficient for a plane wave propagating at the interface between two media for different scenarios. The wave impedance in the air is indicated as Z air , in the skin as Z skin , in the lossless skin as Z lossless skin , and Z effective in the effective medium composed of the nine tissues present in the pectoral model.

		1
		0.9
	Reflection Coefficient | |	0.2 0.3 0.4 0.5 0.6 0.7 0.8		Z effective Z skin ͍ Z air ͍ Z air Z effective ͍ Z lossless skin Z skin ͍ Z lossless skin
		0.1
		0
		10 7	10 8	10 9	5 10 9 10 10
			Frequency (Hz)


Erdem Cil
email: erdem.cil@univ-rennes1.fr

Sema Dumanli
email: sema.dumanli@boun.edu.tr

Denys Nikolayev
email: denys.nikolayev@univ-rennes1.fr

Examination of Impedance Response of Capsule-Integrated Antennas Through Gastrointestinal Tract

Keywords: Capsule antenna, dipole antenna, in-body, ingestible, loop antenna

This work focuses on the investigation of the degree of possible detuning that may be observed in the impedance response of ingestible antennas during the transition through the gastrointestinal tract. For this investigation, two fundamental antennas are used: a meandered dipole antenna and a meandered loop antenna. The antennas conform to the inner surface of the biocompatible capsule shell. They are optimized in spherical homogeneous time-averaged gastrointestinal phantoms to operate at 3 different frequencies of interest: 434 MHz, 1.4 GHz, and 2.4 GHz, resulting in 6 different designs in total. Next, the optimized antennas are simulated in 3 different phantoms, each mimicking the electromagnetic properties of one of the tissues in the gastrointestinal tract (stomach, small intestine, and large intestine). The impedance response of each design in 4 different tissues is compared and discussed. The results show that the maximum shift in the operating frequency of the antennas is 7 MHz, indicating that the presence of the encapsulation layer makes the considered antennas robust against the changes in the electromagnetic properties encountered through the gastrointestinal tract.

I. INTRODUCTION

In recent years, rapid growth has been observed in the use of wireless in-body devices in medical applications. Ingestible devices are one type of these wireless in-body devices. They have a physical structure similar to a capsule and are used by swallowing like conventional pills [START_REF] Kiourti | A review of in-body biotelemetry devices: Implantables, ingestibles, and injectables[END_REF]. Due to their miniature size and compact structure, they are frequently benefited from for various medical purposes such as imaging in a biological system and monitoring diverse physiological parameters inside the human body [START_REF] Kiourti | A review of in-body biotelemetry devices: Implantables, ingestibles, and injectables[END_REF], [START_REF] Rajagopalan | Wireless medical telemetry characterization for ingestible capsule antenna designs[END_REF].

Ingestible devices require various components to be integrated into their capsule to operate effectively [START_REF] Yuce | Easy-to-swallow wireless telemetry[END_REF]. Although the design of ingestible devices is multidisciplinary, the ingestible antenna has a key role in the quality of the wireless link formed by these devices. Hence, the ingestible antenna must be designed meticulously so that it is able to meet the requirements of the wireless link in terms of power and efficiency. Various types of antenna designs have been proposed in the literature for ingestible systems such as helix antenna [START_REF] Faerber | In vivo characterization of a wireless telemetry module for a capsule endoscopy system utilizing a conformal antenna[END_REF], microstrip antenna [START_REF] Nikolayev | Robust ultraminiature capsule antenna for ingestible and implantable applications[END_REF], coil antenna [START_REF] El Hatmi | A multilayered coil antenna for ingestible capsule: Near-field magnetic induction link[END_REF], and inverted-F antenna [START_REF] Kwon | Design of an antenna for an ingestible capsule endoscope system[END_REF].

When a capsule with an ingestible antenna is swallowed, it is transmitted through the gastrointestinal (GI) tract that consists of various tissues, namely the GI tissues, which have different electromagnetic (EM) properties [START_REF] Gabriel | The dielectric properties of biological tissues: II. Measurements in the frequency range 10 Hz to 20 GHz[END_REF]. Therefore, during the design of an in-body antenna, the EM properties of the surrounding environment must be considered, as they strongly affect the antenna parameters such as the operating frequency and antenna efficiency [START_REF] Bosiljevac | In-body antennas design based on fundamental limits of obtainable power density[END_REF]. As the capsule travels through the GI tract, the varying EM properties of different GI tissues may cause a shift in the operating frequency or also affect the antenna efficiency. Thus, an ingestible antenna designer must take into account the changing surrounding environment to ensure the best wireless link quality.

The question one may ask is which antenna parameters must be prioritized while designing an ingestible antenna. For instance, it can be asked whether the designer should concentrate more on the methods to increase antenna efficiency as in [START_REF] Nikolayev | Optimal radiation of body-implanted capsules[END_REF] or it is more essential to increase the robustness against the changes in the environment. To address this question, the work presented in this paper focuses on the degree of the possible detuning in the operating frequency of ingestible antennas that may result from the varying EM properties of the GI tissues. With this investigation, the paper aims to determine whether robustness against the EM changes in the GI tract must be a priority for ingestible antenna designers. To investigate the degree of detuning, a numerical model that consists of a spherical homogeneous phantom and two different types of antennas (dipole and loop) that are integrated into a capsule is prepared as shown in Fig 1 . The antennas are designed to operate at 434 MHz, 1.4 GHz, and 2.4 GHz inside timeaveraged GI phantoms, resulting in 6 antenna designs in total. These frequencies are selected for the examination since they are among the most used frequencies in biomedical in-body applications [START_REF] Damaj | Implantable antennas for biomedical applications: An overview on alternative antenna design methods and challenges[END_REF]. Next, the optimized antennas are simulated in 3 different phantoms that imitate the EM properties of one of the GI tissues (stomach, small intestine, or large intestine) to observe the degree of the shift in the operating frequency.

This paper is organized as follows. In Section II, the method used for the calculation of the EM properties of time-averaged GI phantoms is explained. In Section III, the antenna models are given. In Section IV, the results are shared and discussed. Finally, the paper concludes in Section V.

II. CALCULATION OF EM PROPERTIES OF TIME-AVERAGED GI PHANTOMS

To investigate the degree of detuning during the transmission through different GI tissues, the operation of the antennas must be first optimized in a phantom with specific EM properties. One phantom that can be used for the optimization is a phantom that has time-averaged EM properties of the GI tract [START_REF] Nikolayev | Antennas for ingestible capsule telemetry[END_REF]. To form such a phantom, the time-averaged EM properties are calculated using the time that the capsule spends in each of the tissues. The transition time values are taken from [START_REF] Rao | Investigation of colonic and whole-gut transit with wireless motility capsule and radiopaque markers in constipation[END_REF]. The EM properties of the GI tissues at 3 frequencies of interest are tabulated in Table I along with the calculated timeaveraged values [START_REF]IT'IS Foundation Dielectric Properties Table[END_REF]. For instance, the time-averaged relative permittivity is calculated as

at 434 MHz following this approach. Note that the oesophagus is not taken into account during the calculation, as the transit time from the oesophagus to the stomach is only a few seconds, which is negligible.

III. THE CAPSULE AND THE ANTENNA MODELS

The investigation of the degree of detuning is performed using two simple antenna models (a dipole antenna and a loop antenna) for 3 different frequencies (434 MHz, 1.4 GHz, and 2.4 GHz), resulting in 6 different designs in total. Each antenna is integrated into a capsule. The capsule is based on the capsule given in [START_REF] Nikolayev | Immune-to-detuning wireless in-body platform for versatile biotelemetry applications[END_REF] and consists of a 50 µm thick polyimide (PI) substrate (ε r = 3.4, δ = 0.002) that is wrapped around a cylinder of epoxy-resin (ε r = 3.1, δ = 0.0037), and a 0.6 mm thick Polyvinyl chloride (PVC) layer (ε r = 2.6, δ = 0.003) that covers the PI substrate and the epoxy-resin as shown in Fig 2a . The dipole and loop antenna models are shown in Fig. 2b and in Fig. 2c, respectively, along with the parametrized dimensions. The antennas are meandered on the outer surface of the PI substrate. Meandering is required to fit the physically long antennas into the capsule. The capsules are placed in the center of a spherical homogeneous timeaveraged GI phantom (explained in Section II) with a radius of 50 mm as visualized in Fig. 1. In each case, the dimensions of the antennas are optimized through numerical analysis using ANSYS High Frequency Structure Simulator (HFSS) [START_REF]ANSYS High Frequency Structure Simulator (HFSS)[END_REF] to adjust the operating frequency to the corresponding frequency of interest. The final values for the parametrized dimensions are tabulated in Table II. Note that N shows the number of turns for the dipole antenna and the number of sub-loops for the loop antenna as shown in Fig.

2 [START_REF] Olaode | Effects of meandering on dipole antenna resonant frequency[END_REF]. It can also be noted that the dipole antennas for 434 MHz and 1.4 GHz are offset-fed dipole antennas. Although a center-fed dipole conventionally exhibits a good impedance matching, the offset-feeding is required in these cases to match the antenna to 50 Ω due to long antenna meandering.

IV. RESULTS AND DISCUSSION

Each antenna optimized in Section III is simulated in 3 phantoms that mimic the EM properties of one of the 3 GI tissues at the corresponding frequency: the stomach, the small intestine, and the large intestine (Table I). These phantoms have the same physical and geometrical properties as the timeaveraged GI phantoms. Fig. 3 shows the simulation results in 3 different GI tissue phantoms and in the time-averaged GI phantoms for all designs. From the results, it can be seen that the antennas operate in their intended frequency bands in any GI tissue. The maximum shift in the resonant frequency from the time-averaged GI phantom case is 7 MHz which is observed in the dipole antenna placed in the large intestine phantom at 434 MHz (straight green line and dashed red line in Fig. 3a). These results can be accounted for by the presence of the PVC layer. The PVC layer prevents the antennas from detuning during the transition through the GI tract as most of the fields created by the antennas are contained in this layer. The results indicate that the robustness of the ingestible antennas against the changes in the surrounding environment can be achieved with a proper choice of encapsulation. Therefore, it can be emphasized that an ingestible antenna designer would rather focus on methods to improve other parameters of the antenna such as efficiency rather than robustness. However, it must be noted that this work uses simplified homogeneous phantoms and does not take into account the air pockets that are likely to exist in a physical GI tract. Although the results obtained in this work provide the designers with an insight into the robustness of the ingestible antennas in presence of the encapsulation layer, they are preliminary results and hence open to further investigation with more realistic models.

V. CONCLUSION

This paper examined the degree of possible detuning that may be encountered in the operation of ingestible antennas during the transition through different tissues in the GI tract. A meandered dipole antenna and a meandered loop antenna that are wrapped around a capsule covered with an encapsulation layer are optimized in spherical homogeneous time-averaged GI phantoms for 3 different frequencies. The results show that the considered ingestible antennas are robust against the changes in the EM properties of the GI tract due to the presence of the encapsulation layer. Our ongoing work aims at the fabrication and characterization of the developed antennas in liquid tissue-equivalent models. The results of these measurements will be presented during the conference.

Fig. 1 .

 1 Fig. 1. The numerical model prepared to investigate the degree of the detuning in the operating frequency through the GI tract. The capsule has a length of 20 mm and a diameter of 8.8 mm and the phantom has a diameter of 100 mm.

Fig. 2 .

 2 Fig. 2. The model of the capsule and the antennas and the parametrized dimensions (a) the model of the capsule (b) the model of the dipole antenna and the parametrized dimensions (c) the model of the loop antenna and the parametrized dimensions.

Fig. 3 .

 3 Fig. 3. The simulation results for the dipole and the loop antenna in 3 different GI tissue phantoms and in the time-averaged GI phantoms (a) results for the dipole at 434 MHz (b) results for the loop at 434 MHz (c) results for the dipole at 1.4 GHz (d) results for the loop at 1.4 GHz (e) results for the dipole at 2.4 GHz (f) results for the loop at 2.4 GHz.

TABLE I EM

 I PROPERTIES OF THE TISSUES IN THE GI TRACT AT 3 FREQUENCIES OF INTEREST AND CALCULATED TIME-AVERAGED VALUES

			Tissue
		Stomach Small Intestine Large Intestine Average
	Transition Time (min)	180	228	1302	*
			f = 434 MHz
	Relative Permittivity	67.2	65.3	62	63
	Conductivity (S/m)	1.01	1.92	0.87	1.02
			f = 1.4 GHz
	Relative Permittivity	63.9	57.1	56.1	57.1
	Conductivity (S/m)	1.44	2.44	1.34	1.5
			f = 2.4 GHz
	Relative Permittivity	62.2	54.5	54	54.9
	Conductivity (S/m)	2.17	3.13	2	2.17

ACKNOWLEDGMENT

This study was supported in part by the Région Bretagne through the project SAD "EM-NEURO" and in part by SATT Ouest Valorisation through the project "E-MOTION".