Laurie E Calvet
email: laurie.calvet@c2n.upsaclay.fr

Paoline Coulson

B Iniguez

K Romanjek

2021 Computer-Aided Design for Thin-Film Transistor Technologies (CAD-TFT)

Keywords: organic electronics, spiking neuron circuits, classification circuits

Spike-based neural networks have the potential to significantly increase the computational power of information technologies compared to other architectures. While many different types of sensors and conventional circuits have been implemented using flexible Thin Organic Large Area Electronics, TOLAE, processes, explorations of artificially spiking neuron circuits are much rarer. Here we numerically investigate neuron circuits in gravure printed scalable OTFT technology and consider how the circuits change with variability Finally, we discuss how the intrinsic variability from processing can be used to advantage for classification tasks. These circuits constitute a promising low cost way to integrate neuromorphic hardware with sensors fabricated using the same methods.

Increases in computing performance have led to significant advances in computing systems, ranging from handheld devices to supercomputers, and have enabled rapid advances in Artificial Intelligence (AI). This progress is closely linked to the continued shrinking of the complementary metal-oxide-semiconductor (CMOS) technology feature size, which has now reached fundamental barriers. To improve capabilities of future AI systems, researchers are exploring new computational paradigms, which can exploit the specific needs of a set of applications and the interactions between hardware, software and algorithms in a computing system. Simultaneous to the interest in approximate computation is the explosion of research in printed electronics, which use either organic or oxide based semiconductor materials, for a growing number real-world applications [START_REF] Chang | A Circuits and Systems Perspective of Organic/Printed Electronics: Review, Challenges, and Contemporary and Emerging Design Approaches[END_REF]. Flexible Thin Organic Large Area Electronics (TOLAE) processes enable lower cost, energy and environmental footprint for sensor systems and can provide excellent compatibility with biological objects through their chemistry and intrinsic flexibility [START_REF] Rivnay | Organic electrochemical transistors[END_REF], [START_REF] Cea | Enhancement-mode ion-based transistor as a comprehensive interface and real-time processing unit for in vivo electrophysiology[END_REF]. In addition to consuming less energy during fabrication, at the end of use the electronic materials on flexible substrates can be easily removed and the substrates can be recycled. Nevertheless, organic technology poses large challenges when used for commercial applications. They are slower and although the technology is starting to mature compared to CMOS, they are much less reliable.

Here we consider how to realize hardware for classification of sensors using organic electronics.

We first consider numerically the implementation of a spiking neuron circuit in organic electronics. Next we consider the impact of variability on its functionality. Finally, we consider how to exploit this variability using a carefully chosen classification.

I. SPIKING NEURON CIRCUIT SIMULATIONS

Our neuron circuit simulations are based on the characteristics of 210 OTFT devices with 19 different width/length dimensions. [START_REF] Fattori | A Gravure-Printed Organic TFT Technology for Active-Matrix Addressing Applications[END_REF] The circuit operates in the deep sub-threshold regime, where device characteristics can be modeled as [START_REF] Kim | A Compact Model for Organic Field-Effect Transistors With Improved Output Asymptotic Behaviors[END_REF]:

𝐼 𝐼 𝑒𝑥𝑝 𝑉 𝑉 1 exp (1)
where S is the subthreshold slope, I 0 and V 0 are respectively the current and voltage when the subthreshold current turns on above the background leakage. Note that while V t could also be used, the use of the turn on current and voltage ensures that the extracted parameters are obtained closer to the circuit operation in the deep threshold. Using an automatic fitting program, we extracted the parameters from all the available devices, where channel widths were fixed at 1 mm and lengths varied from 10 to 600 m. The fitting involved 1) obtaining first a fit to the leakage current, 2) identifyig I 0 and V 0 in the leakage current as the points just prior to the large increase in curent, and 3) fitting to equation ( 1) between I 0 and I = 5 x 10 -9 A, the expected maximum of our circuit, to obtain the subthreshold swing S. An example of the fitting is shown in Fig. 1.

Our goal was to realize a biophysical artificial neuron circuit that could realize many of the behaviours exhibited by biological neurons. We modified a Morris-Lecar inspired neuron [START_REF] Calvet | Spiking sensory neurons for analyzing electrophysiological data[END_REF] to accommodate only pOTFT transistors, as shown in Fig. 2. The model neuron gives rise to a spike train due to the opening and closing of 'K' and 'Na' ion channels and their relative timeing. Using Python and the parameters obtained from the real devices, we then implemented circuit simulations. As a starting point for exploring spiking neurons in organic electronics, we used the mean values from the real devices, and then explored how the variation in these parameters effects the circuit performance.

P R P R V GNa V GK V ss V dd V GL C m C k V m P K
Fig. 1 Example of transistor characteristics and fitting procedure to obtain parameters for the compact model in equation 1.

II. RESULTS

Fig. 3 shows an example of the spike trains resulting from the numerical simulations. This circuit functions differently from the CMOS version [START_REF] Calvet | Spiking sensory neurons for analyzing electrophysiological data[END_REF]. Without an excitation current I exc , the P Na transistor is turned on and the P L and P K transistors are in the off state. As I exc increases, the membrane capacitance C m is charged but there are just small changes in V m . Once I exc is large enough to switch the inverter, however, V m exhibits a sharp rise, the P Na transistor will turn off and the P K will turn on. C m will discharge through P L , which will return V m back to its resting potential and the cycle will start again. However, at large I exc , the leak transistor, will remain on and the action of the P Na transistor will be reduced, and thus the spike trains will be lost. The neuron circuit in Fig. 2 thus has the interesting property of exhibiting a relatively small spike range, that can be used to great advantage for sensor applications to limit the impact of larger fluctuations that might occurs for instance to unwanted noise.

We examined the effect of variability by running a series of Monte Carlo simulations, where the parameters of the compact model for each transistor varied up to a percentage. This allows us to estimate how the spike range of the circuit would change for a given variation in the processing. An example of these results in shown in Fig. 4, where random variations of 20% in these parameters still allow for robustly spiking circuits despite large changes in the output characteristics.

III. CONCLUSIONS AND PERSPECTIVES

We have shown how OTFTs can be used to realize artificial spiking neuron circuits. Nevertheless, the circuits are slow and are likely to exhibit relatively large variabilities in their behaviour. The speed of the circuits would enable them to be useful for classification circuits in a wide array of sensors that do not require high speed circuits. In previous work we have shown how spiking neurons can be used as a readout circuit for EEG sensors with significantly reduced data transmission [START_REF] Calvet | Spiking sensory neurons for analyzing electrophysiological data[END_REF]. In addition, they are perfectly compatible with sensors realized in organic technologies [START_REF] Cea | Enhancement-mode ion-based transistor as a comprehensive interface and real-time processing unit for in vivo electrophysiology[END_REF]. Here we argue in addition, that the variability of the circuits can be used to great advantage for realizing classification circuits.

One possibility is to use the spiking neuron circuits in the context of an Extreme Learning Machine (ELM) classification [START_REF] Huang | Extreme learning machine: Theory and applications[END_REF]. In such an implementation, sensor data would be fed into an array of spiking neuron circuits where each one would include a variation in its parameters due to the processing, similar to what is observed in the simulations in Fig. 4. The output of the neurons could be transformed into a non-spiking values via rate encoding, were the number of spikes for a given feature are 'counted'. Subsequently, these values can serve as the input for an offline classification using a machine learning algorithm or the traditional ELM pseudoinverse solution. (In practice such a classification could be carried out on a small device such as a smart phone.) The natural variability of the fabrication results essentially in the original data being projected onto a higher dimensional space so that the machine learning classification can be much more accurate, as in the ELM algorithm. Note that such an implementation does not necessarily result in an increase in the quantity output data relative to the raw sensor data because the spiking neuron can significantly reduce the quantity of the initial sensor data [START_REF] Calvet | Spiking sensory neurons for analyzing electrophysiological data[END_REF]. One could have in fact the best of two worlds: a projection onto a higher dimensional space and at the same time a reduction in the total amount of data, all the while using processes that are low cost and intrinsically variable.

Fig. 2

 2 Fig.2Circuit diagram of a Morris-Lecar inspired spiking neuron circuit simulated with Python using parameters from organic OTFT transistors. Note that the P Na and the P K transistors are meant to be analogous to the 'Na' and the 'K' ion channels in a Morris-Lecar model. The excitatory current is the input and the output is the membrane voltage V m .

Fig. 3

 3 Fig.3Output of the simulated circuit from Fig.2. At small I exc , the response of the circuit is too weak to exhibit spiking behavior, but and increase in I exc results in robust spiking. Finally, at large I exc the spike trains are suppressed.

Fig. 4

 4 Fig. 4 Monte Carlo simulations showing the spike rate versus the Iexc when the parameters of the transistors were randomly varied by 20%. Just one of the circuits was non-functional, while the others always exhibit spiking behavior, albeit with different ranges in Iexc.


Jonathan Cottet
email: jonathan.cottet@epfl.ch

Olivier Fabrègue

Charles Berger

François Buret

Philippe Renaud

Marie Frénéa-Robin

Olivier Fabregue

MyDEP: a new computational tool for dielectric modeling of particles and cells

Keywords: Dielectrophoresis, Cell trapping, Microfluidics, Bottom-up assembly

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

INTRODUCTION

The term "dielectrophoresis" was first introduced by Pohl in 1951 [START_REF] Pohl | The Motion and Precipitation of Suspensoids in Divergent Electric Fields[END_REF] to describe the motion of dielectric particles due to interaction with a non-uniform electric field. Depending on the frequency of the applied field and on the dielectric properties of the particle and its immersion medium, different polarization mechanisms come into play, the main mechanisms being related to the formation of a double electric layer at the particle/liquid interface and to charge accumulation at interfaces between media of different electrical properties (Maxwell-Wagner interfacial polarization effect) [START_REF] Pethig | Dielectrophoresis: Theory, methodology and biological applications[END_REF]. The force resulting from the interaction between the induced dipole moment and the field gradient is expressed by:

= -∇ = . ∇ (1)
where refers to the electric potential energy and E to the electric field.

For a spherical particle of radius , the induced dipolar moment is given by:

= 4 (2)
where is the Clausius-Mossotti factor:

= * - * * + 2 *
(3) * and * refer to the complex permittivities of particle and medium, which depend on their respective electrical conductivities and relative permittivities and on the field angular frequency :

!
where ! is the relative permittivity, the vacuum permittivity, # ! the electrical conductivity and = 2 with f the frequency.

Development of Eq. 1 leads to the expression of the generalized time-averaged DEP force

(3):

= 2 $%&' (∇) *+, - + ./' ( ) -∇0 + ) 1 -∇0 1 + ) 2 -∇0 2 3
(5

)
where 0 , 0 1 and 0 2 are the phase shifts of the field components in the Cartesian coordinates.

Conventional Dielectrophoresis

If the electric field applied is stationary, Eq. 5 simplifies to:

4 = 2 %&' (∇) *+, - (6)
This phenomenon is sometimes referred to as "conventional dielectrophoresis", abbreviated cDEP. The force depends on the gradient of the squared electric field intensity and exists only in the presence of a non-uniform electric field. It is proportional to the volume of the particle, as well as to the real part of the Clausius-Mossotti factor, %&' (. This term, reflecting the polarizability contrast between the particle and its immersion medium, also determines the direction of the force:

• When the particle is more polarizable than its immersion medium %&' ( > 0 , the force acts in the gradient direction and therefore drives the particle towards areas of maximum field intensity. This corresponds to positive dielectrophoresis (pDEP).

• On the contrary, when the particle is less polarizable than its immersion medium %&' ( < 0 , the force moves the particle against the gradient, towards the regions of minimum field intensity, which is referred to as negative dielectrophoresis (nDEP).

Electrorotation

While conventional dielectrophoresis is based on the use of stationary electric fields, electrorotation, abbreviated ROT, induces a rotary motion on a polarizable particle exposed to a rotating field. This technique was developed in the 1980s by Arnold and Zimmermann (4), who proposed to use a four-pole electrode structure to generate the field by applying 90° phaseshifted signals between two adjacent electrodes.

When a polarizable particle is suspended in a rotating electric field, a dipole forms and should rotate synchronously with the field. In practice, when the angular frequency of the field is sufficiently high, the dipolar relaxation time is too long to allow this synchronism. The temporal shift (or phase delay) between the dipole and the field results in a torque exerted on the particle, of expression:

⟨9 : ⟩ = × = -4 ./' () -= > (7)
where = > represents the unit vector normal to the electrode plane and the Clausius-Mossotti factor (cf. Eq. 3).

As the particle rotates, it experiences a resistive viscous torque from the surrounding fluid of expression (5):

⟨9 ? ⟩ = -8 A Ω = > ( 8
)
where η is the dynamic viscosity of the medium and Ω is the constant angular velocity of the particle.

At the equilibrium between induced torque and resistive viscous torque, the rotation rate for a spherical particle is given by:

C = - D D 2A ./' () - (9)
The minus sign indicates that the particle rotates against the field direction for ./' ( > 0. Otherwise the direction of rotation is with the field. Curve fitting procedures may be used to obtain the dielectric parameters of a cell, by minimizing the deviation between the experimental ROT spectrum (plot of the rotation rate with respect to the field frequency) and the theoretical spectrum predicted by the appropriate multi-shell model.

Travelling-Wave Dielectrophoresis (TWD)

The Travelling-Wave Dielectrophoresis force, abbreviated TWD, acts on a particle subjected to a travelling electric field. It is related to the phase non-uniformity of the electric field and arises from the interaction of the travelling field with the phase-lagging component of the induced dipole moment. Such a field can be produced by planar electrodes arranged in rows and driven by a polyphase AC voltage. TWD is therefore an analogue of ROT, Eq 5 remains the same, but with electrodes arranged in line, rather than in a circle. The resulting translational force propels the particle along the electrodes, with or against the field direction, depending on whether ./' ( is negative or positive, respectively. In practice, cDEP and TWD effects can be observed simultaneously: while the particle translates, it is either pushed above the electrodes (nDEP) or attracted onto them (pDEP), depending on the sign of %&' (

.

MyDEP Software

Before performing experiments in the lab with DEP, it is useful to predict the particle and cell responses to the electric field. This requires knowledge of particles or cells properties, which can be obtained from the literature, and implementation of equations related to the particle model.

MyDEP is a computational software, programmed in Java, aiming to study dielectrophoretic behavior of particles and cells in a suspended medium. More precisely the software can calculate and display the Clausius-Mossotti factor (real and imaginary parts) used in DEP for different conditions (medium, frequency range, model of particle). It can also calculate the equivalent permittivities and conductivities of particles alone and in suspension in a medium thanks to the Maxwell-Garnett and Hanai equations. Graphs representing crossover frequencies vs electrical conductivity of the medium are also available.

Cell modeling

Most particles, and especially biological cells, are not homogeneous. It is therefore mandatory to model the different layers that constitute them (cell membrane and cytoplasm in particular). Calculating the Clausius-Mossotti factor requires to successively calculate the equivalent permittivities of the inner layers to obtain a homogeneous equivalent particle. The different models implemented in MyDEP are: "homogeneous particle", "single-shell", "twoshell", "three-shell" and "four-shell" are presented in FIGURE 1.

For a cell modeled with a "single-shell" model composed of a cytoplasm surrounded by a cell membrane, the equivalent complex permittivity (7) is:

D E * = D F * G -Hℎ F J + 2 K D F * -D F * D F * + 2D F * L G -Hℎ F J -K D F * -D F * D F * + 2D F * L (10)
The formulation of the complex permittivity for ellipsoids can be found in [START_REF] Kakutani | Electrorotation of non-spherical cells: Theory for ellipsoidal cells with an arbitrary number of shells[END_REF].

Cell suspension

In the presence of particles, the effective permittivity of the suspension D ! * depends on the volume fraction 0 occupied by the particles. It is given by the Maxwell-Garnett mixing equation if the volume fraction 0 < 0.1

D ! * -D * D ! * + 2D * = 0 D * -D * D * + 2D * (11)
which is equivalent, according to [START_REF] Sihvola | Effective permittivity of dielectric mixtures[END_REF], to the direct formulation:

D ! * = D * N O 1 + 30 D * -D * D * + 2D * 1 -0 D * -D * D * + 2D * Q R ( 12
)
or by the Hanai equation [START_REF] Hanai | Theory of the dielectric dispersion due to the interfacial polarization and its application to emulsions[END_REF], theoretically up to 0 < 0.8 [START_REF] Hanai | Electrical properties of emulsions[END_REF][START_REF] Hanai | Analysis of Dielectric Relaxations of W/O Emulsions in the Light of Theories of Interfacial Polarization[END_REF]:

S D ! * -D * D * -D * T S D * D * T U/ = 1 -0 (13)
Unfortunately, there is no direct expression for the calculation of D ! * and the value can therefore be obtained either by solving the cubic equation or by numerical integration with the difference equation of the Hanai equation. The later has been implemented in MyDEP, with the possibility to choose the number of increments. The methodology can be found in [START_REF] Hanai | Dielectric theory of concentrated suspension of shell-spheres in particular reference to the analysis of biological cell suspensions[END_REF][START_REF] Irimajiri | Dielectric modeling of biological cells. Models and algorithm[END_REF].

Software for dielectric modeling

In 1991, Irimajiri and al. ( 14) published a software for dielectric modeling of particles in suspension called "IMPEDANCE ANALYSIS mini", programmed in BASIC. Unfortunately, this program is not available anymore.

Dielectric modeling of particles behavior is often done in MATLAB (2, 3) and some authors have proposed MATLAB-based programs that are available [START_REF] Erdem | A numerical approach for dielectrophoretic characterization and separation of human hematopoietic cells[END_REF]. Those programs are however limited in terms of functionalities implemented, platform dependent and may lack stability through the different releases of MATLAB and require an access to MATLAB. For those reasons a multiplatform, executable and user-friendly program is needed. This software is intended to be used by both the dielectrophoresis community and for teaching activities. It does not require any prior knowledge of the dielectrophoresis equations.

MATERIALS AND METHODS

The software MyDEP is written in Java using the swing API and is freely available as a standalone .jar file for Windows, Mac or Linux at http://doi.org/10.5281/zenodo.1321928. The installation of Java (also known as the Java Runtime Environment or JRE) is required (https://java.com/en/download). The software version number as well as the recommended java version are available on the website https://mydepsoftware.github.io, under the "Getting Started" menu (item "Installation and system requirements").

A local database is embedded within the application using the SQL database engine SQLite.

The Java Database Connectivity (JDBC) API was used to interact with SQLite.

The static website hosted on GitHub Pages was build using Jekyll. A link to the website is available in the software under the "Help" menu (item "Online resources").

conductivity, # = 7.8 /Y//. The DEP force, which is proportional to Re[CM(f)], is always higher in magnitude for the viable yeast cells.

Parameter sweep

All the parameters of the medium and of the different models can be swept linearly or logarithmically between two values. FIGURE 4 illustrates the graph generated by a logarithmic sweep on ten values of # for Jurkat cells [START_REF] Reichle | A new microsystem for automated electrorotation measurements using laser tweezers[END_REF]. The more the electrical conductivity increases, the lower the initial values of the %&' ( and the shorter the frequency range where Jurkat cells experience pDEP. As # increases, the frequencies of the maximum and minimum of ./' ( are shifted to higher values up to the point where ./' ( keeps a positive value.

Cell separation

Depending on their dielectric properties, different cell types can be separated. FIGURE 5 illustrates that HEK cells [START_REF] Zimmermann | A combined patch-clamp and electrorotation study of the[END_REF] and MCF7 (19) have different responses to the electric field with the frequency. The crossover frequency, transition from the nDEP regime to the pDEP regime, are respectively 169 kHz and 65 kHz for the HEK and MCF7 cells in the specified medium.

These two cell populations can be separated based on their electrical properties between these frequencies. In particular at 100 kHz (vertical blue line in FIGURE 5) HEK cells experience nDEP contrary to MCF7 cells, which experience pDEP.

Conductivity and permittivity of cell and suspension

Crossover frequencies

Crossover frequencies correspond to the frequencies at which %&' ( = 0. It corresponds to the transition from a nDEP regime to a pDEP regime and vice versa. For each electrical conductivity # , this value might differ. FIGURE 9 illustrates the evolution of the crossover frequencies with # for a Jurkat cell. In this figure, the lower crossover, in blue, corresponds to the transition from nDEP to pDEP and the upper crossover frequency, in orange, to the transition from pDEP to nDEP. At approximately # = 0.392 Y//, there is only one crossover frequency point which means that, at a higher conductivity, cells only experience nDEP. Crossover frequencies are commonly used as a discriminatory factor between different cell types

Graph export

All the graph generated in MyDEP are fully editable. The font style and size, the color, the legend content, the curve style and size can all be adjusted directly in the interface. The export menu enables the user to directly generate the displayed graph as an image file with the possibility to tune the size and resolution as well as the file format. A CSV file can be generated if additional data processing is required and not already available in the MyDEP software.

CONCLUSIONS

Figure caption

FIGURE 1 Illustration of the different spherical and ellipsoidal cell and particle models implemented in the interface. All the models "homogeneous sphere", "single-shell", "two-shell" and "three-shell" are illustrated with an example. The implemented "four-shell" model is not illustrated here.

FIGURE 6 ,

 6 FIGURE 6, FIGURE 7 and FIGURE 8 illustrate how cells and medium properties influence

FIGURE 2 FIGURE 3

 23 FIGURE 2 Overview of the database explorer (Search) from MyDEP. A click on the desired element from the literature shows the values associated in the interface.FIGURE 3 Example of %&' ( and ./' ( for viable and non-viable yeast cells suspended in a low conductivity medium (# = 7.8 /Y//, D = 78). Data from (16). The black line represents the baseline at 0.

FIGURE 4

 4 FIGURE 4 Example of the evolution of Re'CM f ( and Im'CM f (with a logarithmic sweep on ten values of σ c from 1 mS/m to 1.6 S/m for a Jurkat cell. Data from[START_REF] Reichle | A new microsystem for automated electrorotation measurements using laser tweezers[END_REF]. The black line represents the baseline at 0.

FIGURE 5

 5 FIGURE 5 Example of the different Re'CM f ( for a HEK cell[START_REF] Zimmermann | A combined patch-clamp and electrorotation study of the[END_REF] and a MCF-7 cell (19) in a medium with σ c = 50 mS/m. The crossover frequencies are respectively 169 kHz and 65kHz. These two cell populations can be separated based on their electrical properties at 100 kHz (vertical blue line).

FIGURE 6 D

 6 FIGURE 6 D E , D and D ! corresponding to respectively the equivalent relative permittivity of a HEK cell, the relative permittivity of the suspension medium and the equivalent relative permittivity of the suspension of HEK cells at a volume fraction 0=0.3 (implemented using the Hanai methodology). σ c = 0.156 S/m.

FIGURE 7 #

 7 FIGURE 7 # E , # and # ! corresponding to respectively the equivalent electrical conductivity of a HEK cell, the electrical conductivity of the suspension medium and the equivalent electrical conductivity of the suspension of HEK cells at a volume fraction 0=0.3 (implemented using the Hanai methodology). σ c = 0.156 S/m.

FIGURE 8

 8 FIGURE 8 | E * |/ , | * |/ and | ! * |/ corresponding to respectively the modulus of the equivalent complex relative permittivity of a HEK cell, the modulus of the complex relative permittivity of the medium and the modulus of the equivalent complex relative permittivity of the HEK cells in suspension in the medium at a volume fraction 0=0.3 (implemented using the Hanai methodology). σ c = 0.156 S/m.

FIGURE 9

 9 FIGURE 9 Evolution of the crossover frequencies for a Jurkat cell for 50 conductivities logarithmically spaced between 1 mS/m and 0.5 S/m.

  voltage-and frequency-dependent membrane capacitance caused by structurally dissimilar lipophilic anions. J. Membr. Biol. 221(2):107-121. 19. Henslee, E. A., M. B. Sano, A. D. Rojas, E. M. Schmelz, and R. V. Davalos. 2011. Selective concentration of human cancer cells using contactless dielectrophoresis. Electrophoresis 32(18):2523-2529.

Acknowledgments

The authors would like to thank the Koji Asami and Marc Castellarnau for the data provided to test the software. The Ampere lab would like to acknowledge support from INSERM (Plan Cancer, Physicancer Program, Dynamo project). The authors also acknowledge the support of the Programme d'Avenir Lyon Saint-Etienne (PALSE mobility grant), the labex iMUST (ANR-10-LABX-0064/ANR-11-IDEX-0007) from University of Lyon as well as the Doctoral school 160 EEA of the University of Lyon for the mobility grant allocated.

RESULTS AND DISCUSSION

Database

MyDEP allows the user to specify the electrical as well as the geometrical parameters of the investigated particle. No prior knowledge of the equations behind is required to use the software, which makes it interesting for users non-familiar with DEP. A database compiling information from the literature is provided to help the user to start with already existing data. The user can also enrich the database with new information. A local database is embedded within the application using the SQL database engine SQLite. The Java Database Connectivity (JDBC) API was used to interact with SQLite. The provided database contains for each set of data the name of the model, the authors, the title of the article, the journal where it was published, the year of publication as well as the DOI or URL to help the users to identify where the model they are using is coming from. An example of the database explorer is displayed in FIGURE 2.

CM factor

MyDEP allows the user to display different graphs linked to the dielectric properties of particles and cells. Real and imaginary parts of the Clausius-Mossotti factor can be displayed in the interface. Users can get values about the displayed curves in the "Results" panel. In particular the values of the crossover frequencies, corresponding to the frequencies at which %&' ( = 0, are directly accessible as well as the frequencies of the minimum and maximum of ./' (. This graph can be used to determine the direction of the DEP force. As shown in


Ali Oulmas
email: ali.oulmas@isir.upmc.fr

Johan E Quispe

Nicolas Andreff

Stéphane Régnier

Comparing swimming performances of flexible and helical magnetic swimmers

HAL is

Comparing swimming performances of flexible and helical magnetic swimmers

Ali Oulmas 1,2 , Johan E. Quispe 1 , Nicolas Andreff 2 and Stéphane Régnier 1

Abstract-Flexible and helical magnetic microswimmers have been well reviewed in the literature because they could be exploited for envisaged applications such as targeted drug delivery, material removal, and micromanipulation. In this article, scaled-up versions of those robots are introduced to study in detail their maneuverability and dexterity while swimming. The robots were immersed in pure glycerol, thus, reproducing a low Reynolds scenario. The proposed robots were previously optimized, achieving their best performances. These magnetic swimmers have different actuation mode and geometric shape. The experiments assess the performances of these two kinds of robots in terms of rapidity, and steering error while following 3D trajectories in environments with high viscous variations.

I. INTRODUCTION

Nowadays, microrobots are good candidates to accomplish envisioned tasks such as minimally invasive diagnosis, targeted therapies, and cell/microobjects sorting [START_REF] Li | Micro/nanorobots for biomedicine: Delivery, surgery, sensing, and detoxification[END_REF]- [START_REF] Bradley J Nelson | Microrobots for minimally invasive medicine[END_REF], [START_REF] Quispe | Exploiting the transient behavior of thermocapillary convection flows to enhance non-contact mesoscale manipulation[END_REF], [START_REF] Quispe | Single particle manipulation/sorting through the transient response of thermocapillary convection flows[END_REF], [START_REF] Zhao | Artificial micro-cinderella based on self-propelled micromagnets for the active separation of paramagnetic particles[END_REF]. To achieve those tasks, a robot suitable to the environment anisotropy and the flow changes must be manufactured. In the biomedical context, these tasks or missions will be in biofluids through human conduits. Moreover, those fluids can pose additional constraints to the robot mission because of their fibrous networks [START_REF] Ullrich | Swimming characteristics of helical microrobots in fibrous environments[END_REF]. Hence, it is preferred to look for a robot capable to overcome every issues presented. Notwithstanding, the above considerations, there exists biofluids such as the cerebrospinal fluid (CSF) that have a viscosity value closes to the water viscosity one [START_REF] Ig Bloomfield | Effects of proteins, blood cells and glucose on the viscosity of cerebrospinal fluid[END_REF] where could serve swimming techniques available for a large interval of low Reynolds number.

As consequence, this work aims at giving a comparison of the capabilities of the helical and flexible magnetic robots at milliscale by reproducing the low Reynolds condition using pure glycerol, and constrained the fabrication in such a way to obtain a free buoyancy millimeter robot. Both robots are tested in various settings considering their widespread control laws, namely giving a homogeneous rotating magnetic field for the helical swimmer, or an oscillating magnetic field for the flexible one. The settings are defined by different trajectories in different viscous environments, and the efficiency is measured in terms of rapidity to execute the trajectories, accuracy to accomplish properly the paths as nicolas.andreff@femto-st.fr well as the manufacture constraints. These criteria are important for microobjects deliverance in order to accurately avoid obstacles and residues in those environments; Analysing the performance of the robot in viscosity changes can help to choose a swimmer which is able to work in a wide variety of (bio-) fluids with different viscosities; Finally, manufacturing process is determinant when it is required to change the swimmer features in order to optimize or adapt its shape for a specific task.

To enfranchise the challenges related to the microscale such as the microfabrication and the perception, a solution consists in studying macroswimmers at low Reynolds number configuration. To do so, it is possible to use a fluid with high viscosity such as glycerol or silicon oil. In this way, standard optical devices such as cameras can be used for 3D observation and tracking and the swimmers are easier to manufacture and customize. These swimmers serve as test benches before validating at microscopic scale. The experimental comparison results confirm the intuition that helical swimmers are more suitable to perform straight line trajectories than flexible ones. Conversely, the flexible swimmers are suitable to perform complex curves with wide range of fluid viscosities.

In the outline of this paper, section II presents the experimental set-up and the different magnetic actuation modes of the helical and flexible swimmers. In section II as well, kinematics equations involved in the system and the control law for the 3D motion control are briefly recalled. Sections III and IV are the core of this work. They define respectively the manufacturing and performances criteria and exhibit the results of helical and flexible swimmers comparison. Finally, section V gives conclusions and new perspectives for this work.

II. EXPERIMENTAL SET-UP AND CONTROL

The experimental system consists of a nested Helmholtz coils system capable of 10mT to actuate the swimmer robots. These robots are provided of neodymium permanent disc magnets in their heads (cf. Fig. 2). Following the procedure in [START_REF] Abbott | How should microrobots swim?[END_REF], Navier-Stokes equation can be written in a dimensionless fashion in such a manner that the flow pattern only depends upon the Reynolds number. In [START_REF] Purcell | Life at low Reynolds number[END_REF], Purcell studied the limit condition when Reynolds number approaches zero. At this limit, inertial forces are negligible compared to viscosity forces, and microorganisms must perform a nonreciprocal pattern in order to swim. The authors in [START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF] reviewed this topic and gave an idea of the magnitude order of the Reynolds number of such microorganisms. For instance, in water, E. coli bacteria achieve a Re of about 10 -5 to 10 -4 , a human spermatozoon a Re of approximately 10 -2 , the Paramecium a Re of around 10 -1 , and so on. In this article, the Reynolds number achieved by both robots is between 10 -3 and 10 -2 . Since these millimetric robots are at low Reynolds, the comparison with their real counterparts at microscale is feasible. With respect to the magnetic actuation, as magnetic field produced by the coils is practically uniform, the robots undergo only a magnetic torque T m as follows:

T m = M × B ( 1
)
where M is the magnetic moment and B the external magnetic field. The latter can be decomposed into driving and steering magnetic fields (respectively B d and B s ) as follows:

B = B d + B s (2)
The driving and steering magnetic field expressions will be detailed in the next sections for both swimmers. However, it is necessary to note that for characterization part only the driven magnetic field was necessary, being for helical propulsion ||B d || = 8 mT, and ||B d || max = 10 mT for flexible propulsion. But also can be written in this way, for helical robots B 0 = 8 mT, and for flexible ones B x = 2 mT and B y = 8 mT, as will be seen.

A. Helical swimmer magnetic actuation

The helical swimmer is actuated by a rotating magnetic field as shown in Fig. 1. The driving and steering magnetic fields are respectively perpendicular and parallel to the principal axis of the robot. The magnetic torque tends to align the magnetic momentum of the swimmer with the applied magnetic field. Thus, with a rotating magnetic field B d and a helical tail, the swimmer propels itself converting its rotation into a linear movement thanks to the force generated by the interaction with the viscous fluid. The expression of the driving magnetic field B d for the helical robot is given as in [START_REF] Oulmas | Closed-loop 3d path following of scaled-up helical microswimmers[END_REF] by:

B d = B 0 cos(2πf t)y B + B 0 sin(2πf t)z B ( 3
)
where f is the rotation frequency, B 0 is magnetic flux density in the center of the work-space while y B and z B are the base vectors of the plane perpendicular to the swimmer axis. The robot steering is made by B s , which is defined by a proportional control of the geodesic error between the real axis of the swimmer and the desired axis [START_REF] Xu | Characterization of three-dimensional steering for helical swimmers[END_REF]: where λ is the proportional control gain and x * B the desired orientation of the helical robot.

B s = -sign(B d • x * B ) λ x B × x * B x B (4)

B. Flagella motion-based actuation

The flexible swimmer is actuated with an oscillating magnetic field. The tail twists when it follows the magnetic field and generates a propulsive force through the interaction with the viscous fluid. The propulsion principle is shown in Fig. 1. The driving magnetic field expression B d is given as in [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF]:

B d = B x x B + B y cos(2πf t)y B ( 5
)
where x B and y B define the swimming plane. Due to the fact that the linear movement direction is not changing with the oscillating magnetic field, a static and homogeneous magnetic field B x = B x x B is added in order to control the flagella orientation. A sinusoidal magnetic field B y = B y cos(2πf t)y B is applied in the direction perpendicular to the axis of the flexible robot, where f is the oscillation frequency. The resulting magnetic field B d oscillates around the flexible swimmer axis. To achieve a desired orientation in space, the steering magnetic field B s is developed as follows:

B s = λ M x * B ( 6
)
where λ is the proportional control gain. The swimmers motion control in closed-loop is briefly introduced in II-C.

C. Modeling and control

The swimmer robots are modeled using a body-fixed frame F B = {x B y B z B } located at the robot center of mass p G with x B the principal axis of the robot. The orientation of the swimmer with respect to the global frame is characterized by the direction angle θ d (yaw) and the inclination angle θ i (pitch). The swimmer principal axis is related to these angles as follows:

x B = [Cθ i Cθ d Cθ i Sθ d Sθ i ] (7)
Fig. 3: Block diagram of the 3D path following for any magnetic swimmer robot behaving as a nonholonomic system. In order to change the propulsion mode specific to the swimmer, only the magnetic controller block is adapted.

This work uses an algorithm based on the path following problem, the aim is to minimize the distance and orientation errors between the swimmer robot and the reference path. Thereby, the state vector is defined as (s, d y , d z , θ de , θ ie ) where s is the curvilinear abscissa of p S (the orthogonal projection of the swimmer onto the path), d y and d z are the projection of the distance between the swimmer and the path on two axes perpendicular to the path tangent vector while θ de and θ ie are the two angles defining the orientation error between the swimmer velocity and the path tangent vector. The kinematics of the swimmers robot in the Serret-Frenet frame is given in [START_REF] Oulmas | Chained formulation of 3d path following for nonholonomic autonomous robots in a serret-frenet frame[END_REF]- [START_REF] Oulmas | 3d closed-loop motion control of swimmer with flexible flagella at low reynolds numbers[END_REF]. To design the controller, an appropriate change of coordinates introduced in [START_REF] Oulmas | Chained formulation of 3d path following for nonholonomic autonomous robots in a serret-frenet frame[END_REF] allows us to transform the swimmer kinematic model into a 5states/3-inputs cascade system from which we were able to derive the following state feedback controller:

u 2 = -k d1 ṡz 2 -k t1 ṡ z 3 (8a) u 3 = -k d2 ṡz 4 -k t2 ṡ z 5 (8b)
where k •• are positive control gains and z •• the state variables.

The swimmers are actually driven by the steering angular velocity Ω = (0 Ω y Ω z ) T , which can be computed as a function of the control inputs:

Ω z = γ -1 21 (u 2 -γ 22 ) (9a
)

Ω y = γ -1 31 (u 3 -γ 33 -γ 32 Ω z ) (9b)
where γ •• are appropriate scalars depending on the states and the path [START_REF] Oulmas | Chained formulation of 3d path following for nonholonomic autonomous robots in a serret-frenet frame[END_REF]. Finally, the desired orientation x * B can be computed thanks to the time integration of the steering angular velocity:

x

* B = Ω × x B dt ( 10
)
where dt is the sampling time. To summarize, the block diagram of the entire algorithm, used to control the robot motions in space, is presented in Fig. 3. The path following algorithm (i.e. kinematic controller block on the figure) uses only the spatial errors between the robot and the reference trajectory, independently of the robot type. Thus, it is sufficient to adapt the actuation mode specific to the swimmer (i.e. the magnetic controller block) to switch from the helical swimmer motion control to the flexible one.

III. COMPARISON OF TOPOLOGICALLY DIFFERENT

SWIMMERS

Helical and flexible swimmers use both non-reciprocal movements for displacing inside liquid media at low Reynolds number. However, they are completely different geometrically. Thereby, processes involved in the manufacturing are an important difference to consider. In this work, the helical swimmer has a rigid helical-shaped tail, fabricated using a 3D printer with "VisiJet M3 Black" material. Moreover, a disc-shaped magnet with 1.5 mm in diameter and 0.5 mm in height was attached to the swimmer head in such a way that the magnetization is perpendicular to the swimmer principal axis [START_REF] Oulmas | Closed-loop 3d path following of scaled-up helical microswimmers[END_REF]. The second swimmer has a flexible tail, made from an elastomer (silicone) by mixing base and catalyst liquids, which cures at ambient temperature [START_REF] Oulmas | 3d closed-loop motion control of swimmer with flexible flagella at low reynolds numbers[END_REF]. A certain ratio between the flexible swimmer elasticity and the fluid viscosity has to be respected. This relation is described by the sperm number:

S p =
Viscosity forces Elastic forces [START_REF] Oulmas | Closed-loop 3d path following of scaled-up helical microswimmers[END_REF] For S p >> 1, the tail is very flexible and its oscillations decrease rapidly because of the high viscosity value of the fluid. In the other case, S p << 1, the tail is rather rigid and it does not twist during the oscillations. In this situation, the swimmer movement is reciprocal and the "scallop theorem" is applied, thus, the swimmer velocity tends to zero. For the magnetization, the same disc-shaped magnet is attached to the flexible tail in such a manner to have the magnetization in the same direction as the robot principal axis. Some characteristics of the actual helical and flexible swimmers are summarized in TABLE I. Concerning the manufacturing, helical swimmers are straightforward to fabricate at millimeter scale rather than flexible swimmer because they can be manufactured using a 3D printer which allows to make complex-shaped and various swimmers. On the other hand, flexible swimmer requires to accomplish with the flexibility ratio or sperm number (Sp) that means to take into account the material properties such as its Young modulus, its Poisson ratio, among others. In our particular case, we made use of a mold where base and catalyst are mixed strictly in 1:1 proportion in order to respect the ratio between the tail elasticity and the fluid viscosity. Besides, currently the magnetization method (by surface functionalization or by providing a magnet) for both robots has not been completely mastered, especially at microscale. As the manufacturing of the helical swimmer is mastered at milliscale, we propose an experimental optimization of the helical swimmer geometry in III-A. Different swimmer shapes are proposed.

A. Helical swimmer optimization

Helical swimmers have been well reviewed in the literature [START_REF] Gao | Bioinspired helical microswimmers based on vascular plants[END_REF]. Among the most relevant features, the axial pitch, length, number of turns and magnet shape are studied in various works [START_REF] Xu | Characterization of three-dimensional steering for helical swimmers[END_REF]. However, the number of tails has not been considered yet in the case of rigid tail. In this subsection is presented an optimization with several helical tails considering a fixed length. In order to know how the number of tails affects the swimmer dynamics, different helical swimmers are designed (cf. Fig. 4) by varying the number of tails and mass. Fig. 4 shows that the forward speed of helical swimmers depends on the rotation frequency and the number of tails. It is well known that the helical robot reaches its maximal value of forward speed at the cut-off or stepout frequency [START_REF] Gao | Bioinspired helical microswimmers based on vascular plants[END_REF], [START_REF] Xu | The rotational propulsion characteristics of scaled-up helical microswimmers with different heads and magnetic positioning[END_REF]. From our new results, under the conditions described in the caption of Fig. 4, we can observe that the swimmer speed decreases by increasing the number of flagella. This could be caused by the hydrodynamic interactions of the flow surrounding the different flagella when the space between them becomes smaller. Consequently, the total speed of the swimmer decreases. Thus, in this paper, helical swimmers with single flagellum present better performances in terms of propulsion and will be used in the rest of comparison experiments below. Conversely to the intuition, this result show that by increasing the number of flagella, speed decreases. Different to results demonstrated for multiple flexible flagella [START_REF] Zhou Ye | Rotating magnetic miniature swimming robots with multiple flexible flagella[END_REF].

IV. COMPARISON OF SWIMMERS WITH DIFFERENT PROPULSION MODES

In this work, two kinds of bio-inspired swimmers are designed, manufactured and their movements are controlled in closed-loop. However, which swimmer to choose to make specific trajectories in particular environments? This section tends to answer this question by comparing the swimmers experimentally using the curvature radii and the fluid viscosity.

A. Radius of curvature criterion

The capacity of the swimmers to perform complex trajectories can be evaluated by studying their performances in terms of accuracy and execution rapidity to different curvature radii of the curves to be realized. In actual applications, this criterion can help to choose the swimmer which is able to perform complex paths in environments full of obstacles and residues (such as blood) by keeping high accuracy. In this work, the helical and flexible swimmers are nonholonomic, which means that they are limited for performing paths with complex curvatures. For that, a series of experiments with circular trajectories of different radii are performed to evaluate their abilities. The paths realized by the flexible swimmer are depicted in Fig. 5 with curvature radii ranging from 1 mm to 5 mm with a pitch of 1 mm. For the helical swimmer, the curvature radii of the circles are varying from 2 mm to 10 mm with an axial pitch of 2 mm instead of 1 mm as for the flexible swimmer in order to maintain the same ratio length/radius for each swimmer. In fact, the helical swimmer length is twice as large as the flexible swimmer length. The paths performed by the helical robot are shown in Fig. 5.a. To study the precision of these two swimmers, the root mean square errors (RMS) of the lateral distance d y during the path following of different circles are reported in Figure 5.d. For the flagella, the error is less than 500 µm in all cases. For the helix, the error is 1 mm for a radius of curvature of 2 mm and for the other curvature radii, the accuracy is quite satisfactory. Thus, if we take the swimmers at the actual size, the results show that the flexible swimmer is much more precise than the helical one. The performances of these swimmers are also studied in terms of execution rapidity. TABLE II shows the time that the two robots take to perform a circle turn. It can be seen that except for a normalized radio of 12.5, the helical swimmer is faster than the flexible one because it takes less time to make a turn of circle. To conclude, based on the experiments of this section, helical swimmers are more convenient to perform tasks with straight line trajectories because they are faster than the flexible ones but less accurate when the trajectory curvature becomes complex. On the other hand, flexible swimmers are more suitable for paths with complex curvature because they are accurate but very slow comparing to the helical swimmers. Furthermore, as the helical swimmer is twice as long as the flexible swimmer, a comparison can be made by normalizing the results and the curvature radii of the reference paths by the length of the swimmers, as illustrated in Fig. 5.c. The graph shows that for the same length, the precision performance of the helical swimmer could be better than the flexible one. To summarize, at the current state with the manufacturing processes that we master, the performance of the flexible swimmer is better than the helical robot in terms of accuracy during the path following of complex curvatures.

B. Fluid viscosity criterion

The swimmers dynamics depend strongly on the environment in which they are evolving, particularly the fluid viscosity. Their cut-off (or optimal) frequencies can change in function of the viscosity since the fluidic torque depends on it. To study the ability of the swimmers to operate in fluids with different viscosities, experiments have been performed using vertical sinusoidal trajectories in closedloop with two different viscosities of glycerol solution. First, a pure glycerol (99.99 % with viscosity of 1.5 Pa s) is used which is referenced as a high viscosity fluid for swimmers. Then, a 1% in volume of water is added, reducing the fluid viscosity to 20 % from the initial value (from 1.5 to 1.2 Pa s [START_REF] Sheely | Glycerol viscosity tables[END_REF]). This is referred as a low viscosity fluid in these experiments. The rotation and oscillation frequencies of the magnetic field are kept fixed to the maximal values for all the experiments, as well as, the control gains.

For high viscosity, the results are depicted in Figure 6.a. It can be seen that both swimmers converge towards the reference path and follow well the trajectory. The distance quadratic errors of both swimmers once on the path (X ≈ 23mm) are less than 200 µm. For comparison, the same are performed with a lower fluid viscosity. In actual applications, the fluid is not always homogeneous. Therefore, it is important to evaluate the behaviours of these swimmers in different conditions such as different fluid viscosities. The results of the path following are illustrated in Figure 6.b. It is noticed that the flexible swimmer converges towards the desired trajectory even in low fluid viscosity. However, it can be clearly seen that the helical swimmer has difficulties to reach the desired path. In these conditions, the vertical velocities for both swimmers are measured for different inclination angles and are shown in figure 6.c. We can see that for the flagella, the vertical velocity increases proportionally with the inclination angle. However, the helix vertical velocity reaches a plateau after 30 • at a weak value. This explains why the helical robot can not follow the desired trajectory. From these experiments, we can conclude that the flexible swimmer would be potentially able to work in wide range of fluids with different viscosities comparing to the helical swimmer. These results are valid for the robots used in this paper. After swimmers optimization, these results may be different, however, the fluid viscosity remains an interesting criterion of comparison between robots with different topologies. Fig. 7 gathers different criteria discussed in this paper for the flexible and helical swimmers comparison. These criteria are tested experimentally and help to choose the most efficient swimmer for a given trajectory and viscosity.

V. CONCLUSION This paper proposes an experimental comparison between the bio-inspired helical and flexible swimmers on complex trajectories using multiple criteria such as the curvature radii of the desired trajectories and the fluid viscosity. The presented results allow to choose the most efficient swimmer for a given fluid viscosity and reference trajectory. Thus, helical swimmers are more suitable to perform straight line trajectories and the flexible ones for complex curves with wide range of fluid viscosities. This study attempts to propose a set of criteria for comparing swimmers that present different topologies and propulsion modes. The aim is to provide guidelines for users to choose the most efficient robot relative to a specific task in the future. As the fluid motion around the microswimmer is characterized only by the Reynolds number, the propulsion at the microscopic scale should be similar to that of the macroscopic swimmers. As a next step, microswimmers will be investigated to validate this work at low scales. For that, the macroswimmers will be used as guidelines, for example, the ratios length/diameter and magnetic torque/fluidic torque will be kept the same to get the same performances. In addition, these comparison criteria will be used similarly to choose the most performant microswimmer as this study is dimensionless. As perspectives, it would be interesting to exploit more criteria such as boundary effects which can affect the swimmers motion. This measure can be important for tasks that take place in a confined space. The consumed energy criterion can be also studied.

Fig. 1 :

 1 Fig. 1: Flexible and helical swimmer and their actuation mode by oscillating and rotation uniform magnetic fields in (a) and (b) respectively.

Fig. 2 :

 2 Fig. 2: The experimental set-up for the wireless magnetic manipulation of the helical and flexible swimmers at low Reynolds numbers.

Fig. 4 :

 4 Fig. 4: On the right side, forward speeds vs rotation frequency in pure glycerol for the multi-tail helical-shaped robots. On the left side the manufactured robot and its respective CAD model for 1, 2, 3 and 4 tails in (a), (b), (c) and (d) respectively.

Fig. 5 :

 5 Fig. 5: Radii of curvature criterion. Circle tracking of the flexible (a) and helical (b) for different radii ranging from 1 to 5 mm for the flexible and from 2 to 10mm for helical, constant gains were kept constant in both cases. (c) Lateral error graph normalized by the swimmer length. (d) Lateral distance RMS errors of the path following of the helical and the flexible.

Fig. 6 :

 6 Fig. 6: Viscosity criterion. (a) Sinusoidal tracking trajectory of helical and flexible swimmer in the vertical plane ZX using a high viscosity fluid and results for low viscosity are depicted (b). (c) Depicts the vertical velocity of the helical and flexible swimmers using the low viscosity fluid.

Fig. 7 :

 7 Fig. 7: Helical and flexible swimmers comparison with different criteria.

TABLE II :

 II Time spent by the swimmers to complete a turn in seconds. The curvature radii are normalized by the swimmer length. Hel and Fla present respectively the abbreviations of the helical and flexible swimmers.

	Normalized radii	12.5 %	25 %	37.5 %	50 %	62.5 %
	Time (s)	Hel. Fla.	54.01 34.1	78.15 85.51	104.73 132.12	135.06 141.23	148.46 203.71

*This work was supported by EIPHI Graduate School (ANR-17-EURE-0002), MultiFlag (ANR-16-CE33-0019-02) and Grand Prix Scientifique 2018 Ch. Defforey-Institut de France.