Valerie Cognat

Gael Pawlak

David Pflieger

Laurence Drouard
email: laurence.drouard@ibmp-cnrs.unistra.fr

Valérie Cognat

PlantRNA 2.0: an updated database dedicated to tRNAs of photosynthetic eukaryotes

Keywords: transfer RNA genes, pol III transcription, plant genomes

INTRODUCTION

Beyond their major role in translation as carriers of amino acids, transfer RNAs (tRNAs) have many non-canonical functions and are involved in a multitude of biological processes, as already emphasized in PlantRNA 1.0 database [START_REF] Cognat | PlantRNA, a database for tRNAs of photosynthetic eukaryotes[END_REF] or in more recent reviews (e.g. [START_REF] Su | Noncanonical Roles of tRNAs: tRNA Fragments and Beyond[END_REF] [START_REF] Tosar | Editorial: Understanding the Importance of Non-Canonical tRNA Function[END_REF]). For instance, the key role played by tRNAs in the production of tRNA-derived small RNAs (tDRs) [START_REF] Holmes | A Standardized Ontology for Naming tRNA-derived RNAs Based on Molecular Origin[END_REF] involved in diverse functions is expanding rapidly [START_REF] Krishna | tRNA-derived fragments (tRFs): establishing their turf in post-transcriptional gene regulation[END_REF]. Furthermore, there is growing evidence that not only tRNAs but also tRNA genes appear to play important roles. In particular, tRNA genes can be involved in the organization of the nuclear genome and affect the chromatin structure (e.g. [START_REF] Van Bortle | tDNA insulators and the emerging role of TFIIIC in genome organization[END_REF] [START_REF] Hamdani | tRNA Genes Affect Chromosome Structure and Function via Local Effects[END_REF] [START_REF] Hummel | The multi-faceted regulation of nuclear tRNA gene transcription[END_REF]) [START_REF] Ehrlich | On the Track of the Missing tRNA Genes: A Source of Non-Canonical Functions?[END_REF]). These genes are either scattered or clustered throughout the genomes [START_REF] Bermudez-Santana | Genomic organization of eukaryotic tRNAs[END_REF] [START_REF] Hummel | Epigenetic silencing of clustered tRNA genes in Arabidopsis[END_REF]. Their number greatly varies between the nuclear genomes of eukaryotes [START_REF] Hummel | The multi-faceted regulation of nuclear tRNA gene transcription[END_REF], and likely their genomic distribution varies as well but this remains poorly explored. To better decipher how tRNA genes contribute to genome evolution it is therefore essential to retrieve complete sets of tRNA genes with the highest possible accuracy.

Photosynthetic organisms, with the presence of three genomes (i.e. nuclear, plastidial, and mitochondrial) represent the most complex model of eukaryotic cells to visualize the evolution of the tRNA genes population in terms of gene content, organization and expression. Different tRNA gene databases including plant tRNA genes have been developed over time (Table 1). We can cite the Transfer RNA database (tRNAdb), [START_REF] Juhling | tRNAdb 2009: compilation of tRNA sequences and tRNA genes[END_REF], a restructured version of the first compilation of tRNA and tRNA gene sequences created in 1989 [START_REF] Sprinzl | Compilation of tRNA sequences and sequences of tRNA genes[END_REF], the tRNA Gene DataBase Curated by Experts (tRNADB-CE), [START_REF] Abe | tRNADB-CE 2011: tRNA gene database curated manually by experts[END_REF], the Genomic tRNA Database (GtRNAdb), [START_REF] Chan | GtRNAdb 2.0: an expanded database of transfer RNA genes identified in complete and draft genomes[END_REF], and the T-psi-C database (T-Psi-C) [START_REF] Sajek | T-psi-C: user friendly database of tRNA sequences and structures[END_REF]. In 2011, we developed PlantRNA, the first tRNA database dedicated to photosynthetic organisms and provided information on 11 eukaryotes of the photosynthetic lineage [START_REF] Cognat | PlantRNA, a database for tRNAs of photosynthetic eukaryotes[END_REF]. Although covering evolutionary distinct branches of the photosynthetic lineage, this database was rather limited. During the last decade, the number of completed nuclear, plastidial and mitochondrial genome sequences has considerably increased. This recent accumulation of knowledge has now made it possible to include the tRNA gene sequences of about 40 additional species covering the phylogenetic tree of the Archeplastida lineage (PlantRNA 2.0, http://plantrna.ibmp.cnrs.fr/). This includes species representative of terrestrialisation, vascularisation, flowers emergence, or growth in extreme environments as well as two rhodophytes whose genomes contain permuted tRNA genes or peculiar intronic sites. The genomes of a cryptophyte, a haptophyte and two heterokontophytes were also annotated. The recently optimized version of tRNAscan-SE program [START_REF] Chan | tRNAscan-SE 2.0: improved detection and functional classification of transfer RNA genes[END_REF] was used to predict and annotate tRNA genes as accurately as possible. Nevertheless, manual curation was still essential at different levels, for example, to avoid false annotations due to contamination from external (e.g. bacterial) or organellar DNA during nuclear genome sequencing programs, to eliminate remaining tRNA pseudogenes that escape detection, or to identify missing plastidial tRNA genes containing group II introns. Biological information relevant to tRNA biology (e.g. intron sequences, flanking sequences controlling tRNA gene expression, tRNA mitochondrial import) is also made available. To note, two other tRNA databases dedicated to plants, PltRNAdb (https://bioinformatics.um6p.ma/PltRNAdb/index.php ), [START_REF] Mokhtar | PltRNAdb: Plant transfer RNA database[END_REF] and PtRNAdb (http://14.139.61.8/PtRNAdb/index.php) [START_REF] Singh | PtRNAdb: a web resource of plant tRNA genes from a wide range of plant species[END_REF] were published in 2022, demonstrating the need to compile plant tRNA gene sequences.

As a whole, circa 37 000 novel tRNA genes were accurately annotated, thus providing the largest highly curated tRNA database of photosynthetic organisms.

RESULTS

Database content and web interface

The current release incorporates data sources of the previous version for 7 species [START_REF] Cognat | PlantRNA, a database for tRNAs of photosynthetic eukaryotes[END_REF], novel data sources for 4 updated genomes, and a choice of 40 new data sources selected on the basis of the interest of the species regarding their position in the phylogenetic tree of the photosynthetic lineage. For instance, they cover 7 genomes of streptophyte algae belonging to the lineage being at the transition step from aquatic algae to land plants, the genomes of a lycophyte (Selaginella moellendorffii) the first non-seed vascular plant genome reported (Banks et al., 2011), a fern (Azolla filiculoides), one of the fastest-growing plants, or Amborella trichopoda, representative of the sister lineage of all flowering plants. Plant species were also chosen because of their growth in peculiar environments such as the cave plant Primulina huajiensis [START_REF] Feng | The genome of a cave plant, Primulina huaijiensis, provides insights into adaptation to limestone karst habitats[END_REF] or the halophytic plant Eutrema parvulum (also called Schrenkiella parvula or Thellungiella parvula) [START_REF] Dassanayake | The genome of the extremophile crucifer Thellungiella parvula[END_REF][START_REF] Wang | The Whole Genome Assembly and Comparative Genomic Research of Thellungiella parvula (Extremophile Crucifer) Mitochondrion[END_REF]. The quality of the genome annotations was also taken into account whenever possible. In addition, when available, plastidial and mitochondrial genomes were also used as sources. Whole nuclear, plastidial, and mitochondrial genomes were scanned by tRNAscan-SE 2.0 [START_REF] Chan | tRNAscan-SE 2.0: improved detection and functional classification of transfer RNA genes[END_REF] and automatically processed by a homemade workflow (available at https://github.com/ibmp-bip/tRNAflow) and then manually curated and annotated as in [START_REF] Michaud | A global picture of tRNA genes in plant genomes[END_REF]. Figure 1 shows the phylogenetic tree of the photosynthetic organisms from which the tRNA genes were retrieved. The sources of genome sequences listed in Supplemental Table S1 are accessible via the entry point "Species" on the PlantRNA website. As in version 1.0 [START_REF] Cognat | PlantRNA, a database for tRNAs of photosynthetic eukaryotes[END_REF], additional data including linear secondary structure, A and B boxes, flanking sequences controlling tRNA gene expression, intronic sequences, mitochondrial tRNA import either experimentally proven or predicted, are indicated. Regarding intronic sequences, introns predicted at non-canonical positions in the red alga, Chondrus crispus, were experimentally verified. To this end, total RNA was extracted from C. crispus alga provided by the 'Centre national de ressources biologiques marine' (EMBRC France) according to (Salinas-Giege et al., 2017) and analyzed by RT-PCR or cRT-PCR (circular RT-PCR) as described in [START_REF] Warren | Rapid Shifts in Mitochondrial tRNA Import in a Plant Lineage with Extensive Mitochondrial tRNA Gene Loss[END_REF] and using oligonucleotides found in Supplementary Table S2. We also annotated some of the sequences as tRNA-like sequences. This tRNA-like group comprises true tRNA pseudogenes but likely also highly divergent tRNA molecules which functions remain to be elucidated. Without experimental studies, it is difficult to discriminate between the two categories. In addition, numerous mitochondrial or plastidial tRNA gene sequences inserted into nuclear genomes, and reported so far not be expressed and functional [START_REF] Michaud | A global picture of tRNA genes in plant genomes[END_REF], or contaminating the nuclear DNA preparations are recognized as true tRNA genes by tRNAscan-SE. They were annotated and can be downloaded in xls format. Finally, tRNA genes of bacterial origin likely due to DNA contamination were also discarded. All tRNA sequences, annotations, and biological information are stored in a database implemented in PostgreSQL version 12 (https://www.postgresql.org/).

To increase the traceability, we added the genome resources (genome assembly version and source) and the software and versions used for automatic annotation. The web application querying the underlying database is running in a Java-based Apache Tomcat 8.5 environment (https://tomcat.apache.org/). To manage the increasing number of genomes to process and improve the data quality, some automation and systematic data checks are included in a private area of the application.

Global picture of tRNA gene populations

In total, 34403, 1661, and 911 tRNA genes were registered from the nuclear, plastidial, and mitochondrial genomes. Among the 51 photosynthetic organisms, 4 plastidial and 13 mitochondrial genome sequences are not available. At the level of nuclear genomes analysis, the level of curation was variable (Figure 2). tRNA genes of bacterial origin retrieved on scaffolds (i.e. not localized in chromosomes) by tRNAscan-SE 2.0 were deleted: 47 in Coffea canefora, 40 in Mesotaenium endlicherianum, 35 in Mesostigma viride, and 20 in Chondrus crispus. Similarly, many plastidial or mitochondrial tRNA genes were detected; their number greatly fluctuates from 0 in the hornwort Anthoceros angustus to 1978 in the fern Azolla filiculoides. When located on chromosomes, the presence of these tRNA genes is due to organellar DNA insertion into the nuclear genomes; when unplaced (i.e. on non-assembled scaffolds), it is difficult to discriminate between organellar DNA insertion and organellar DNA contamination, in particular when genomes have not been assembled into chromosomes. Finally, depending on the quality of the genome assembly and annotation, the number of tRNA-like/tRNA pseudogenes is highly variable. More than 10 000 sequences from the nuclear genomes of the fern A. filiculoides or of the gymnosperm Pinus taeda belong to this category as well as circa 6 000 sequences in the genome of the streptophyte alga, Chara braunii. As shown in Figure 2, the remaining true tRNA genes may represent not more than 10% (e.g. in P. taeda or in C. braunii) of the initial sequences retrieved thanks to tRNAScan-SE. The presence of numerous tRNA-like/tRNA pseudogenes may be explained by genome size expansion for example by whole-genome duplication as is the case in ferns or in gymnosperms for instance [START_REF] Li | Fern genomes elucidate land plant evolution and cyanobacterial symbioses[END_REF] [START_REF] Zimin | An improved assembly of the loblolly pine megagenome using long-read single-molecule sequencing[END_REF]. In gymnosperms, the high proportion of tRNA pseudogenes is in agreement with the numerous pseudogenes observed in conifer genomes and is correlated with the accumulation of transposable elements [START_REF] Nystedt | The Norway spruce genome sequence and conifer genome evolution[END_REF]. All photosynthetic organisms studied here possess intron-containing nuclear tRNA genes. In embryophytes (Figure 1), around 5% of them possess a classical intronic sequence between positions 37 and 38 and correspond to two tRNA gene families (tRNA Tyr and tRNA eMet ). This proportion increases up to 13.8% in Brassicaceae, such as A. thaliana, where a cluster of tRNA Tyr genes exist [START_REF] Hummel | Epigenetic silencing of clustered tRNA genes in Arabidopsis[END_REF]. In gymnosperms, Ginko biloba and Pinus taeda, tRNA Ile (TAT) genes also contain classical intronic sequences. In streptophyte algae and chlrorophytes, the number and identity of intron-containing nuclear tRNA genes is more variable and usually higher. The permuted tRNA genes identified in the rhodophyte Cyanidioschyzon merolae [START_REF] Soma | Permuted tRNA genes expressed via a circular RNA intermediate in Cyanidioschyzon merolae[END_REF] were also annotated as well as the intronic sequences located at non-canonical positions in the D-, T-or variable loops [START_REF] Kawach | Unique tRNA introns of an enslaved algal cell[END_REF]. Similarly, scanning the nuclear genome of another rhodophyte, Chondrus crispus, several essential tRNA genes were considered as pseudogenes when using tRNAScan-SE. Manual inspection allowed us to identify introning sequences in either the D-loop between positions 20 and 21, or in the T-loop between positions 58 and 59 of these tRNA genes. Interestingly, the introns in the D-or T-loops are always at the same position and tRNA Tyr and tRNA eMet genes have two introns, the canonical one in the anticodon region plus a second one in either the D-or T-loop (Figure 3). The expression of these C. crispus tRNAs and the splicing of their introns were supported by RT-PCR and cRT-PCR experiments conducted for three of them, tRNA Cys (with an intron in the D-loop), tRNA Trp (with an intron in the T-loop), and tRNA Tyr (with two introns in the D-and anticodon loops). In the three cases, correctly spliced tRNAs ending with a CCA triplet were obtained (Supplementary Figure S1). Regarding mitochondria-encoded tRNAs, the novel annotated genomes extend our knowledge on the number and identity of this population of tRNA molecules. Interestingly, even only considering embryophytes, there is a great fluctuation in their number: from 0 in the lycophyte Selaginella moellendorffii to 126 in the basal angiosperm Amborella trichopoda.

To compensate for the lack of mitochondria-encoded tRNAs, this likely implies mitochondrial import of a full set of nuclear-encoded tRNAs in the lycophyte, a process experimentally demonstrated in several photosynthetic organisms [START_REF] Warren | Rapid Shifts in Mitochondrial tRNA Import in a Plant Lineage with Extensive Mitochondrial tRNA Gene Loss[END_REF]. Conversely, the high number of mitochondria-encoded tRNAs in A. trichopoda reflects the presence of numerous horizontal transfers of foreign mitochondrial DNA from green algae, mosses and angiosperms [START_REF] Rice | Horizontal transfer of entire genomes via mitochondrial fusion in the angiosperm Amborella[END_REF]. Even with the existence of 126 mitochondria-encoded tRNA genes, the only tRNA Thr (GGT) is likely not sufficient to decode the four threonine ACN codons, thus the need to import at least one nuclear-encoded tRNA Thr is expected. Based on the lack of mitochondriaencoded tRNA species, imported tRNAs were proposed in the PlantRNA database for all sequenced mitochondrial genomes.

DISCUSSION

Nuclear tRNA genes are transcribed by RNA polymerase III (pol III). The high number of available nuclear tRNA sequences now allow studying in more detail the important characteristics of pol III transcription. In particular, it will be possible to follow how conserved are the TATA-like elements and CAA motifs found upstream of tRNA gene sequences in photosynthetic organisms and whether their presence is linked to plant evolution or adaptation to the environment. Similarly, stretches of T residues normally involved in pol III transcription termination are not always present or are variable in length or number. Whether this variability reflects differences in genome organization and expression remains to be elucidated. As mentioned in the introduction, six other databases, where tRNA genes from photosynthetic organisms have been annotated, are freely available (Table 1). The oldest one, tRNAdb only contains 247 annotated plant nuclear tRNA gene sequences, mostly from A. thaliana, and no full genome has been analyzed. This database has not been updated since 2009. Nonetheless, it is an essential database, as it provides not only tRNA gene sequences but also tRNA sequences with important information on nucleotide modifications, including 55 plant tRNA sequences. More recently, the T-Psi-C database has been released. It includes an updated collection of tRNA sequences obtained through high-throughput tRNA sequencing and many 3D tRNA structures. This tRNA database is particularly relevant for scientists working on human cells for instance. However, only very few plant tRNA sequences are available, mainly those found in tRNAdb. In the third tRNA database, tRNADB-CE, the identified tRNA genes were manually checked by experts and thus provide highly reliable sets of tRNA genes. A huge number of genomes were analyzed (in particular, numerous metagenome analyses of microorganisms), including 2325 complete plastidial genomes but only two nuclear genomes of photosynthetic organisms were annotated and the database has not been updated since 2011. The fourth genomic tRNA database, GtRNAdb, is a very nice and valuable follow-up of the powerful tRNAscan-SE prediction program on complete or nearly complete genomes and 17 land plant genomes were analyzed. Thanks to the efficiency of the tRNAscan-SE 2.0 tool, the data provided are very useful regarding angiosperms studies and score of confidence are easily accessible. However, no lower photosynthetic organisms are integrated in the analysis and no organellar genomes annotated. Very recently, two tRNA gene databases dedicated to photosynthetic organisms were made available: PtRNAdb and PltRNAdb. Using tRNAscan-SE 2.0 and various filters to report reliable tRNA sequences, the PtRNAdb compiles tRNA gene sequences from 106 nuclear genomes mainly from angiosperms (all except 3), and from 89 and 38 plastidial and mitochondrial genomes respectively. However, due to the absence of a manual curation taking into account published data, some information is missing or incorrect. Non exhaustively, we can cite the lack of discrimination between initiator and elongator tRNAs, false annotation of organellar tRNA Ile (CAT) as tRNA Met [START_REF] Weber | A potato mitochondrial isoleucine tRNA is coded for by a mitochondrial gene possessing a methionine anticodon[END_REF], or the absence of nucleotide G-1 present in organellar tRNA His genes and in cytosolic tRNAs His [START_REF] Rosen | Evolutionary conservation of a functionally important backbone phosphate group critical for aminoacylation of histidine tRNAs[END_REF] [START_REF] Placido | Plant mitochondria use two pathways for the biogenesis of tRNAHis[END_REF]. Quite importantly, some intron-containing plastidial tRNA genes are missing (because not retrieved by tRNAscan-SE), while additional mitochondrial tRNA pseudogenes are wrongly annotated. As an example with the model plant A. thaliana, 4 plastidial tRNA genes (Gly, Val, Leu and Lys) are omitted [START_REF] Sato | Complete structure of the chloroplast genome of Arabidopsis thaliana[END_REF] and its mitochondrial genome only encodes 22 tRNA genes [START_REF] Unseld | The mitochondrial genome of Arabidopsis thaliana contains 57 genes in 366,924 nucleotides[END_REF] not 35 as retrieved in PtRNAdb. Thus, although valuable and very useful, data provided in PtRNAdb have to be carefully checked. In the PltRNAdb database, nuclear and organellar tRNA gene sequences from 256 plant species were identified. This database is very extensive with plants belonging to Streptophyta, Chlorophyta, and Rhodophyta. For tRNA genes prediction, tRNAscan-SE 2.0 was used but no other curation has been performed and relevant published literature has not been taken into account. Taking A. thaliana as a case study, as above for ptRNAdb, plastidial intron-containing tRNA genes are lacking, the number of 36 mitochondrial tRNA genes exceeds by far the true number of 22, organellar tRNA Ile (CAT) is annotated as tRNA Met gene and the G-1 of tRNA His sequences is missing. Such types of errors exist for all plant species. As other examples, circa 130 tRNA Leu genes retrieved from the Populus trichocarpa genome are not true tRNAs Leu , the introns present in C. crispus tRNAs were not predicted, the Chlamydomonas reinhardtii tRNA Sec [START_REF] Rao | Chlamydomonas reinhardtii selenocysteine tRNA[END_REF] is missing and the set of tRNA genes of bacterial origin contaminating the Coffea canephora genome has not been deleted. As a whole, ptRNAdb, represents an extensive database for a rapid overview of the tRNA gene populations found in photosynthetic organisms. Here, the PlantRNA database provides tRNA gene sequences from nuclear and organellar genomes of only 51 photosynthetic organisms but these organisms were selected to reflect the whole phylogenetic tree and our stringent manual curation ensures the recovery of tRNA gene sets of high quality where relevant biological information have been considered. Thanks to these two criteria, choice of species and high curation, the PlantRNA database is complementary to the other databases and represent an excellent resource for researchers in tRNA biology. In PlantRNA 2.0, no emphasis has been given on related biological information such as aminoacyl-tRNA synthetases or modification enzymes and work needs to be done to fill this gap. From an informatics point of view, many steps required to implement the PlantRNA database have now been automatized. The main bottleneck in updating more regularly the database remains the essential manual curation. We now plan, with the large number of curated tRNA sequences, to implement a neural network so that any new genome is accurately annotated and the data quickly incorporated into the database. Finally, the incorporation into the database of tRNA and tDR expression profiles as well as the modified nucleotide composition of each tRNA species is biologically relevant but will first require more experimental work.

. Phylogeny of Archaeplastida with the species included in this study. Two heterokontophytes belonging to the SAR (Stramenopiles, Alveolates, Rhizaria) clade are included. The tree has been constructed from various published studies, in particular (Amborella Genome, 2013, [START_REF] Burki | The New Tree of Eukaryotes[END_REF][START_REF] Rensing | How Plants Conquered Land[END_REF][START_REF] Wang | Genome-wide analyses across Viridiplantae reveal the origin and diversification of small RNA pathway-related genes[END_REF]. LECA: Last Eukaryotic Common Ancestor. For each species, a three letters abbreviation is indicated in bold.

Figure 2. Histogram showing the percentage of nuclear tRNA genes, as estimated from our study. Species are indicated by the three letters code presented in Figure 1. Note that two Chlamydomonas reinhardtii strains (Cre-503 and Cre-4532) have been studied. The color code is as followed: true tRNA genes in blue, tRNA-like/tRNA pseudogenes in orange, tRNA genes inserted into the nuclear genome or due to contamination by mitochondrial or plastidial DNA in grey, and tRNA genes corresponding to external bacterial contamination in yellow.

Figure 3 .

 3 Figure 3. Insertion points introns in Chondrus crispus tRNA genes. The typical secondary structure of a mature tRNA is presented. The anticodon is represented with an orange circle. Insertion points of canonical introns in the anticodon loop, and of non-canonical introns in either the D or the T-loops are indicated by arrows. The tRNAs possessing introns are listed by their cognate amino acid, their anticodon in parenthesis, and the number of isodecoders.

  1365313x, ja, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/tpj.15997 by Université De Strasbourg, Wiley Online Library on [14/10/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Table 1 . List of tRNA database.

 1 , ja, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/tpj.15997 by Université De Strasbourg, Wiley Online Library on [14/10/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Species 1365313x, ja, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/tpj.15997 by Université De Strasbourg, Wiley Online Library on [14/10/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

	100%
	90%
	80%
	70%
	60%
	50%
	40%
	30%
	20%
	10%
	0%
	A t h	A ly	C r u	B n a	B r a	B o l	E p a	E s a	P t r	P v u	G m a	M t r	M d o	F ii	R c h	V v i	C c a	A m a	P h u	S ly	S t u	C a n	H a n	H v u	B d i O s a Z m a	A t r	P t a	G b i	A f i S m o	P p a M p o	A a n P m a M e n	S m u	C b r	K n i	C a t	M v i C r e -5 0 3 C r e -4 5 3 2	O t a C m e	C c r	C p a	G t h	D lu	E s i P t r i

1365313x

1365313x, ja, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/tpj.15997 by Université De Strasbourg, Wiley Online Library on [14/10/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

ACKNOWLEDGEMENTS

We would like to thank Marie Bernhard and Thalia Salinas-Giegé for technical support, and Alexandre Berr for critical discussion.

FUNDING This work was supported by the Centre National de la Recherche Scientifique (CNRS) and by the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm supported by IdEx Unistra (ANR-10-IDEX-0002), and EUR IMCBio (ANR-17-EURE-0023) under the framework of the French Investments for the Future Program.

1365313x, ja, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/tpj.15997 by Université De Strasbourg, Wiley Online Library on [14/10/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

DATA AVAILABILITY STATEMENT

All data can be found within the manuscript, its supporting materials and via the freely accessible PlantRNA web interface (http://PlantRNA.ibmp.cnrs.fr). All genome resources (genome assembly version and source) and the software and versions used for automatic annotation are indicated in the manuscript and available on the website. All published data performed with the help of the PlantRNA database should refer to this article.

AUTHOR CONTRIBUTIONS

LD conceived the project, curated manually the data, and wrote the manuscript. VC developed the workflow to annotate tRNA, processed all the genomes, and contributed to the design of the database and the data integration. GP implemented the PlantRNA database and the web interface and DP adapted the tRNAscan import R package to split the tRNA sequences. All authors contributed to the final manuscript.

CONFLICT-OF-INTEREST STATEMENT

No potential conflict of interest was reported by the authors.

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article. S1. List of the sources of genome sequences Table S2. Primers used in this work.


Giacomo Fanciullo

Irene Conti

Pascal Didier

Andrey Klymchenko

Jérémie Leonard

Marco Garavelli

Ivan Rivalta

Jérémie Léonard

Modelling quenching mechanisms of disordered molecular systems in the presence of molecular aggregates

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

Introduction

The theoretical modelling of the exciton density dynamics allows to shed light on the photophysical processes happening within photoactive materials, such as those employed in the development of artificial light-harvesting systems and organic light emitting diodes. Typical materials of interest are those based on conjugated polymers in the solid phase [1][2][3][4], where structural defects localize the excited states into small portions of the polymer chain acting as independent dyes, but several studies have been performed on aggregates in the solid phase [5,6], thin solid films of organic dyes [7][8] and on the so-called host-guest systems, where organic dyes are spatially distributed within polymeric hosts in form of films [9,10] or nanoparticles [11,12]. If dyes are weakly interacting, one can assume excitons localized on individual dyes and thus able to move through the dye distribution. This exciton mobility can be described by invoking an incoherent energy transfer (ET) hopping mechanism, which can be successfully described by the Förster theory of the ET [13,14]. In Förster model, the timescale governing the donor-acceptor ET is strictly dictated by the radiative (or fluorescence) lifetime of the donor molecule. Typical lifetimes for fluorophores range within the nanosecond scale, but several ETs have been experimentally observed to happen also on the picosecond [2,15] and ultrafast femtosecond [16] timescales.

Figure 1.

The mechanism for exciton-exction annihilation (EEA). The excitation is transferred between two excited dyes (DYE 1 and 2), leading to an excited state absorption (to the N th excited state) in one dye (DYE 2) and to a relaxation to ground state in the other one (DYE 1). The N th high-lying is expected to fastly decay non-radiatively to the first excited state, leading to the annihilation of one excitation (i.e. that of DYE 1).

Other than the exciton mobility, the ET can be responsible for quenching mechanisms which remove excitations from the system. A quenching process driven by ET is the so-called exciton-exciton annihilation (EEA), which leads to a loss of excitations from the system due to the interaction between the excited states of two molecular dyes in spatial proximity. As shown in Figure 1, the EEA occurs when an excited dye transfers its excitation energy to another excited dye, promoting it to a higher excited state while returning to its ground state. The lifetime of higher excited states is usually very small due fast non-radiative internal conversion processes which are likely to relax back the dye to its first excited state. Since the higher excited state can be rapidly depopulated via non-radiative mechanisms, the whole process describes a fast loss of excitation from the system.

The time evolution of the exciton density, which can also occur in an ultrafast time scale [17], is generally studied by the mean of transient absorption [2,3,[5][6][7][8][9]17] or photoluminescence decays [10,18,19] measurements. The recorded, time-resolved signal is proportional to the time-dependent exciton density within the sample, and the proportionality constant is usually recovered by determining the number of absorbing dyes, i.e. their excitation probability under illumination by the light pulse used. In photoluminescence decays experiments, straightforward evidence for the presence of EEA is the acceleration of the fluorescence signal decay upon increasing the excitation power : since powerful irradiations excite a higher number of dyes, excited dyes become closer in space thereby increasing their annihilation probability. Still, the decay of the excitation density can be also affected by other competing quenching phenomena, such as the ET to nonfluorescent aggregates of dyes that is often neglected, while the role of molecular aggregation in the photophysical properties of optical systems is attracting increasing attention. For instance, the effect of the aggregation has been extensively studied in the application of the rhodamine B based dyes employed in the building of artificial light harvesting nanosystems [11,20]. Rhodamine B aggregates are known to be fluorescence quenchers [21,22] and it has been demonstrated that the use of bulky hydrophobic counterions efficiently prevent aggregation [23,24]. Previous theoretical investigations suggested that such quenching mechanism could be due to the fact that rhodamine B dimers (with H-type aggregation) feature dark charge-transfer excited states very close in energy to the bright states, allowing for internal conversion processes that provide suitable paths for non-radiative decays. [22] Reliable theoretical simulations of the fluorescence decay acceleration in presence of EEA can thus provide essential information on the spatial distribution of molecular excitations and on the concomitant role of molecular aggregation on photoactive materials. The theoretical simulation of the time evolution of the exciton density in presence of EEA could be derived through a microscopic description based on the quantum master equation formalism [25], which gives the time evolution of the population of each individual dye, but this level of theory requires a deep knowledge of the dye distribution, in particular of the distances between dyes and their relative orientation. In fact, such approach has been applied to conjugated polymers for which simplified but reliable spatial models are available [26][27][28], but for disordered systems detailed information about the structure cannot be straightforwardly achieved. A simpler level of description invokes for a macroscopic description of the exciton density by strictly assuming a homogeneous distribution of dyes. Exclusively in this framework, the exciton density 𝑛(𝑡) (quantifying the number of excited dyes per volume unit) in presence of EEA can be described using the following equation:

𝑑𝑛(𝑡) 𝑑𝑡 = - 1 𝜏 𝑛(𝑡) - 1 2 𝑘 ,,-(𝑡) 𝑛(𝑡) / (1)
Here, the first term describes the decay due to the excited state dye lifetime, with lifetime 𝜏 comprising both radiative and non-radiative processes, while the second describes the decay due to EEA. As shown in Figure 2, three different models have been developed to express the rate constant 𝑘 ,,-. Please do not adjust margins Please do not adjust margins

In the "direct mechanism" model, the EEA is assumed to happen as described by the Förster ET theory (Figure 2A), and the rate constant takes the form originally derived by Forster [13] for a generic donor-acceptor pair:

𝑘 ,,-(𝑡) = 2 3 𝜋𝑅 1 * 1 * 7 8 𝜋 𝜏 9:1 𝑡 (2)
where 𝑅 1 * 1 * is the Förster radius related to the ET between two excited dyes and 𝜏 9:1 is the fluorescence lifetime of the dye. The 𝑡 ;<// time dependence indicates the fact that excitations in close proximity annihilate firstly, thus reducing with time the probability of having annihilation events. In the "diffusive mechanism" model, the excitations are allowed to diffuse through dyes, and the EEA is assumed to take place when two excitations reach a certain contact radius 𝑅 4 (figure 2b). The rate constant takes the form derived for the generic problem of a diffusing particle captured by an immobile capturing center [29]:

𝑘 ,,-(𝑡) = 8𝜋𝐷𝑅 4 ?1 + 𝑅 4 √2𝜋𝐷𝑡 B (3)
Here, 𝐷 is the diffusion coefficient describing the excitation mobility within the dye network. If we consider excitations moving in the dye network through a Förster hopping mechanism, the diffusion coefficient may be expressed as [30][31][32][33]]

𝐷 = 𝜂 𝑅 1 * 1 D 𝜏 9:1 ? 4 3 𝜋𝐶 1 B G/7 (4)
Where 𝑅 1 * 1 is the Förster radius related to the ET between an excited and a relaxed dye, 𝐶 1 is the dye concentration and 𝜂 = 0.43 is a factor accounting for the homogeneous distribution of dyes. The capability of eq.( 4) to correctly estimate the diffusion coefficient has been analysed in detail by Colby et al [10], showing that its successful application could be compromised by several conditions such as the dye aggregation, the breakdown of the dipole-dipole interaction to describe the electronic coupling or the anomalous diffusion driven by energetic and orientational disorder. In this work, we will assume a framework in which the expression for the diffusion coefficient given by eq.( 4) is valid and the dye aggregation contribution can be treated separately (as the aggregates, acting as quenchers, will be not directly involved in the diffusion processes). Eq. ( 4), in fact, allows for an interesting parametrization of the diffusive-related phenomena by simply introducing the additional Förster radius 𝑅 1 * H . Gösele et. al. [34][35][36][37] derived equations for the rate constant considering both direct and diffuse mechanisms, as depicted in Figure 2A and 2B, and their results have been re-obtained starting from a more sophisticated model by Jang et al [38]. They distinguish two different regimes, having one among the diffusive or the direct ET mechanisms as predominant. For each regime, they also consider two limiting temporal scales: the short-times scale, where only one mechanism is assumed to rapidly start, and the long-times scale, where a stationary regime is assumed, i.e. the probability of having a decay due to both mechanisms is considered constant in time. The authors proposed that an approximate form for the global rate constants (one for each regime) should be obtained by simply summing the expressions for the two temporal scales. In the predominant diffusion regime, at short times only the diffusive mechanism is assumed to happen, and the global rate equation takes the same form of the pure diffusive case, eq.( 3). In this work we will focus on the predominant direct ET regime (see Figure 2C), which has been already suggested for simulations of the exciton density decays [9]. In this regime, the diffusion is assumed not to be started yet at short times, and the rate constant 𝑘 ,,- IJK9L (𝑡) becomes the same of the Forster model in eq.( 2), while at long times the rate constant takes the form

𝑘 ,,- MKNO = 4𝜋𝐷 0.676 T 𝑅 1 * 1 * D 𝜏 9:1 𝐷 U </G (5)
It is worth noting that although the direct mechanism is predominant, for the rate constant exists a dependence on the diffusion coefficient. The global rate constant is finally calculated as

𝑘 ,,-(𝑡) ≈ 𝑘 ,,- IJK9L (𝑡) + 𝑘 ,,- MKNO (6)
which can be rewritten in a compact form as

𝑘 ,,-(𝑡) ≈ 4𝜋𝐷 𝑅 ,,-? 1 + 𝑅 ,,- √𝜋𝐷𝑡 B ( 7
)
where a new "annihilation radius" 𝑅 ,,-has been introduced

𝑅 ,,-= 0.676 T 𝑅 1 * 1 * D 𝜏 9:1 𝐷 U </G (8)
Since the expression obtained for the rate constant has the same form of the diffusive mechanism one, eq.( 3), we can consider the combination of the two processes as a new "effective diffusion" process, with a new contact distance equal to the annihilation radius (see Figure 2C).

As mentioned above, it could be the case that EEA is not the only process responsible for the decay of the excitation density: if the system contains intruder quenching species (such as non-fluorescent aggregates of dyes), the ET from excited dyes to these molecular quenchers will act as a new channel for wasting excitations. Once the excitation energy has been transferred to quenchers, they likely do not relax instantaneously but rather feature a non-radiative lifetime, which implies that an excited quencher density evolving in time has to be considered. In this work, the excitation density This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins Please do not adjust margins time-evolution of a photoactive system involving both the EEA and the dye-to-quencher ET phenomena (with dyes' aggregates as quenchers) is phenomenologically modelled based on the model of Gösele et al. in a regime where the direct ET is the predominant effect. We report a detailed analysis of the effects of the presence of a small amount of quenchers and of the influence of all model parameters used to simulate the exciton density decays in disordered molecular systems. In order to describe the applicability of the model to a realistic system, we consider as a hypothetical case of study the distribution of alkyl rhodamine B dye with bulky counterions within the organic nanoparticles presented in Ref. [11], from which we take some parameters, such as the dye's fluorescence lifetime and concentration. We finally discuss on the usefulness of model for the global fitting of experimental data and its applicability, which should stimulate future developments to obtain model parameters from firstprinciples simulations.

Results and discussion

The time evolution of the excitation densities for dyes 𝑛 1 (𝑡) and quenchers 𝑛 H (𝑡) (we assume no triplet states population) is described by the following system of coupled equations: where 𝑛 H YZY -𝑛 H (𝑡) is the density of relaxed quenchers expressed in terms of the total concentration of quenchers 𝑛 H YZY , and with 𝜏 1 and 𝜏 H being the total lifetimes of dye and quencher, respectively (accounting for both non-radiative and radiative processes in the case of the dye).

𝑑𝑛 1 (𝑡) 𝑑𝑡 = - 1 𝜏 1 𝑛 1 (𝑡) -𝑘 ,,-(𝑡) 𝑛 1 (𝑡) / -𝑘 W (𝑡)𝑛 1 (𝑡)X𝑛 H YZY -𝑛 H (𝑡)[ (9)
The ET processes considered in our model are depicted in Figure 3. The rate constants for the EEA are described by eq.s ( 7) and ( 8) whit a diffusion coefficient given by eq.( 4), and similar expressions are considered also for the dye-quencher ET:

𝑘 W (𝑡) = 4𝜋𝐷 𝑅 H ? 1 + 𝑅 H √𝜋𝐷𝑡 B ( 10
)
𝑅 W = 0.676 T 𝑅 1 * H D 𝜏 1 9:1 𝐷 U </G (11)
In the following, the analysis of parameters affecting the numerical solution of eq.( 9) is reported by assuming a slight amount of quenching species within the dye distribution. We consider as reference dye rhodamine B, having a radiative lifetime 𝜏 1 9:1 = 4.24 ns (corresponding to the measured value for the 0.1% loaded nanoparticles in ref. [20]) with a concentration of 𝐶 1 YZY = 0.13 dyes/nm 3 (corresponding to the 30% loaded nanoparticles in ref. [20]). We initially assume the quencher not to exceed 1% of the dye concentration, thus simulating a typical experimental situation in which a residual aggregation persists in the system despite the efforts made to prevent it. We neglect, for simplicity, the possibility of having non-radiative contributions to the lifetime of both dye and quencher. Equations ( 9) has to be solved imposing certain initial condition, that is the exciton density at the time zero, 𝑛 1 (0): we consider that immediately after the light pulse that 7% of the total dyes are excited (a population value falling in the typical experimental range [9]), while we assume that no quenchers are initially excited. This assumption is based on the fact that a small amount of aggregates is present in the sample (we generally consider <1% of 𝐶 1 YZY ) and likely their absorption maximum is shifted from that of isolated dyes since, to work as quenchers, a strong electronic coupling has to be present among the monomeric units. The effect of varying the initial values of the exciton density simulates the typical acceleration of the exciton density decay upon power light increase in experiments (an example is provided in the SI, see Figure S1 in section S.1). The following analysis aims highlighting the various effects of the parameters entering our model, emphasizing those features that could allow a straightforward identification of the presence of intruder quenching species. We also explored how much such analysis could depend on the initially chosen values for dye/quencher lifetimes and concentrations, showing the adaptability of our tools to various molecular systems with ET processes occurring at different time scales, from nanoseconds to the ultrafast.

Please do not adjust margins

Please do not adjust margins The influence of the 𝑅 1 * 1 radius is quite different from that of 𝑅 1 * 1 * , as clearly shown in Figure 5. Differently from the radius governing the direct EEA, it affects both the short and long times evolution. While the influence on the short times is similar to 𝑅 1 * 1 * , here also the long-time decays are accelerating upon increasing the radius. We now consider the effect of the parameters related to the quencher species, i.e. its concentration 𝐶 H , Förster radius 𝑅 1 * H and radiative lifetime 𝜏 1 9:1 , as reported in Figures 6, 7 and 8A, respectively. It is evident that when the quencher is present in a small amount, its influence can affect the long times evolution only, since the few quencher molecules are very far from the majority of excited dyes and then a certain time has to pass to allow excitations to diffuse towards them. Since this specific (isolated) effect on the long-time evolution cannot be obtained by varying the parameters related to the EEA only, it represents a straightforward signature for the presence of intruder quenching species. Notably, all the trends described above involving an initial amount of excited dyes (𝑛 1 (0) = 7% of the 𝐶 1 YZY ), i.e. Figures 234567, are preserved when such initial (fixed) parameter is significantly increased (see Figures S2-S7 in section S.2 of the SI, for trends with 𝑛 1 (0) = 25% of the 𝐶 1 YZY ). Please do not adjust margins

Please do not adjust margins

It is important to note that the timescale at which the ET processes take place (which in the Förster theory are defined by the molecules' lifetime) does not affect the outcome of the analysis reported above. For instance, as shown in Figure 8, the effect of the quencher lifetime 𝜏 H for ET dynamics in the picosecond and femtosecond timescales (Figures 8B and8C, respectively) exhibits the same features observed in the nanosecond timescale (Figure 8A). The effects of the other parameters (Forster radii and quencher concentrations discussed in Figures 4567) are also the same at different timescales, as reported in section S.3 of SI. The advantage of including quenching species in cases of disordered molecular systems, such as those containing rhodamine B taken as exemplifying case in this work, relies on the fact that a global fitting of experimental decays by using eq.( 9) could demonstrate and eventually quantify the presence of intruder quenchers due to aggregation phenomena. It is worth mentioning that we performed some preliminary tests on real experimental data, observing that the model is able to well fit fluorescence decays when fixing up to two parameters in the proposed model. Clearly, such ability of the model to fit experimental data has to be ascribed to its very high flexibility, originating from the large number of free parameters fitted, including the 𝑅 1 * 1 and 𝑅 1 * 1 * Förster radii (governing the diffusive and direct mechanism of ET processes among dyes, respectively) and the three parameters related to the quencher, i.e. its lifetime 𝜏 H , concentration 𝐶 H and 𝑅 1 * H Förster radius. Thus, to make applications of our model bringing real physical information, one has to be able to experimentally and/or theoretically determine the highest number of them. The 𝑅 1 * 1 can be straightforwardly obtained from the overlap between the steady state absorption and emission spectra of the dye [14,[START_REF]FRET -Förster Resonance Energy Transfer[END_REF]. The 𝑅 1 * 1 * Förster radius could be obtained from the transition absorption spectroscopy by isolating the excited state absorptions contribution and calculating the overlap with the emission spectrum of the dye [9,[START_REF] Nettels | [END_REF]. This is less straightforward than obtaining 𝑅 1 * 1 since excited state absorptions often overlap with ground state bleaching and stimulated emission signals. In this context, advanced simulations of nonlinear electronic spectroscopy from firstprinciples, such those developed in our group in the recent years [41][42][START_REF] Segarra-Martí | Towards Accurate Simulation of Two-Dimensional Electronic Spectroscopy[END_REF][44][45], would be of great support to apply the presented model to real systems. Still, getting reliable deconvolution of transient absorption spectra into specific (overlapping) signals is a very challenging task also from a theoretical point of view, since signal lineshapes are difficult to simulate quantitatively. [42] Finally, regarding parameters associated to the dye aggregates (working as quenchers in our model), first-principles simulations would also be crucial to estimate the 𝑅 1 * H Förster radius, especially if such aggregates cannot be isolated and their absorption spectra recorded. Using structural models of these aggregates, indeed, absorption and emission spectra could be simulated and their

Förster radius thus estimated. In such a way, one could restrict the set of fitting parameters in our model to the quenchers' concentration and lifetime only, providing fundamental insights into the role of aggregation in the photophysical properties of photo-responsive molecular materials. Please do not adjust margins Please do not adjust margins

Materials and methods

The system of coupled equations ( 9) has been solved numerically by means of an in-house Python code employing the scipi.integrate package, available upon request to the authors. The rhodamine B dye assemblies in polymeric nanoparticles of [20] have been used as exemplifying case for a disordered molecular system. We considered particles with 45 nm diameter and 30% w/w dye loading. Form the average number of rhodamine B molecules per particles reported in [20], we calculated the dye concentration 𝐶 1 YZY . For the total lifetime 𝜏 1 , we considered the value measured for the nanoparticles with the lowest loading (0.1% w/w) because in such a diluted sample a negligible aggregation is expected, and the measured lifetime should be representative for that of isolated dyes. The radiative lifetime 𝜏 1 9:1 has been calculated from the fluorescence quantum yield 𝛷 of the low-loaded nanoparticles through the relation 𝛷 = 𝜏 1 /𝜏 1 9:1 .

Conclusions

The exciton density dynamics in photoactive materials is significantly affected by the possibility of having energy transfer (ET) phenomena and the combined use of timeresolved spectroscopy and theoretical modelling has become a widely used tool to recover fundamental information on the type and on the extent of ET processes. When the photoactive system is composed of a uniform distribution of donors and acceptors, simple macroscopic models can be exploited to describe the exciton density dynamics. Such macroscopic models, involving both direct (Förster-like) and diffusive (with ET described by classical diffusion) ET mechanisms, have been successfully applied to describe the influence of excitonexciton annihilation (EEA) on the exciton density dynamics, but they usually neglect any other quenching phenomena competing with EEA. One of these phenomena that is attracting significant attention is the formation of dye aggregates acting as fluorescence quenchers (e.g. rhodamine B and derivatives).

In this work, we proposed a macroscopic model to simulate exciton density dynamics in disordered molecular systems where both direct (Förster-like) and diffusive (via classical diffusion) ET mechanisms are involved, and allowing both exciton-exciton annihilation (EEA) and ET to molecular aggregates as possible quenching mechanisms. To explore the applicability of the model on a realistic system, we considered as study case a homogeneous system in which rhodamine B dyes are embedded in polymeric nanoparticles. We performed a systematic analysis on the role of various parameters, including the Förster radius 𝑅 1 * 1 that is used to express the dye excitation diffusion coefficient, the Förster radii 𝑅 1 * 1 * and 𝑅 1 * H governing quenching mechanisms via the EEA among dyes and dye-to-aggregate ET, respectively, and the aggregate quencher's concentration and lifetime. The results suggested that while the 𝑅 1 * 1 and 𝑅 1 * 1 * inevitably affect the short-time behaviour of the decay, the parameters related to the quenchers have an influence only on the long-time decays, providing a fingerprint of aggregation effects on exciton density dynamics. Moreover, we showed that the effect of a small amount of quencher exhibits always the same features, independently from the power of the light source experimentally used and, more importantly, independently from the ET timescales, including the case of ultrafast dynamics. The overall outcome indicates that the application of our model to real cases could allow straightforward identification of quenching effects arising from molecular aggregates. This could be achieved by global fitting analysis of experimental data that, however, could be compromised by an extreme flexibility of the model, in case of too many parameters are set free. Thus, we envisioned that experimental determination of some of the model parameters could be supported by first-principles simulations. In particular, when the experimental determination is very challenging, as for parameters of the quenching aggregated species or those derived from deconvolution of time-resolved electronic spectra of the dye, first-principles estimates of parameters could enter in our model, potentially providing unique information on the role of molecular aggregation on the photophysical properties of photoactive materials. This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins Please do not adjust margins

Figure 2 .

 2 Figure 2. Schematic representation of possible EEA mechanisms. A) The direct mechanism follows the Förster theory and the ET is governed by the Förster radius 𝑅 1 * 1 * . B) The diffusive mechanism involves a sequence of energy hops moving from the excited dye through the unexcited ones described by a diffusion coefficient 𝐷, while the EEA is assumed to take place when the exciton enters the contact sphere defined by the contact radius 𝑅 4 . C) The model by Gösele et al. (in the predominant direct mechanism regime as proposed in Ref. 8) describes the combined effect of both direct and diffusive mechanisms as an effective diffusive mechanism driven by the usual diffusion coefficient 𝐷 but governed by a new contact radius 𝑅 ,,-, which in turns depends on both the Förster radius 𝑅 1 * 1 * and 𝐷.

  𝑡) + 𝑘 W (𝑡)𝑛 1 (𝑡)X𝑛 H YZY -𝑛 H (𝑡)[

Figure 3 .

 3 Figure 3. The two possible mechanisms for the exciton loss in presence of both EEA and quenching species (such as dimers in the case of rhodamine B) described by the Gösele model in the predominant direct mechanism regime. The new introduced contact radii 𝑅 ,,-and 𝑅 H depends on the Förster radii 𝑅 1 * 1 * and 𝑅 1 * H respectively, and depend also on the same diffusion coefficient 𝐷.

Figure 4 .

 4 Figure 4. Effect of the 𝑅 1 * 1 * radius governing the direct EEA mechanism on the (normalized) exciton density decay. Fixed parameters are: 𝑅 1 * 1 = 5 nm; 𝑅 1 * H = 2 nm; 𝜏 H = 4 ns; 𝐶 H = 0.25 % of 𝐶 1 YZY ; 𝑛 1 (0) = 7 % of 𝐶 1 YZY .

Figure 5 .

 5 Figure 5. Effect of the 𝑅 1 * 1 radius governing both EEA and quenching diffusive mechanisms on the (normalized) exciton density decay. Fixed parameters are: 𝑅 1 * 1 * = 4 nm; 𝑅 1 * H = 2 nm; 𝜏 H = 4 ns; 𝐶 H = 0.25 % of 𝐶 1 YZY ; 𝑛 1 (0) = 7 % of 𝐶 1 YZY .

Figure 6 .

 6 Figure 6. Effect of the 𝑅 1 * H radius governing the direct quenching mechanism on the (normalized) exciton density decay. Fixed parameters are: 𝑅 1 * 1 * = 4 nm; 𝑅 1 * 1 = 5 nm; 𝜏 H = 4 ns; 𝐶 H = 0.25 % of 𝐶 1 YZY ; 𝑛 1 (0) = 7 % of 𝐶 1 YZY .

Figure 7 .

 7 Figure 7. Effect of the quencher concentration 𝐶 H (expressed as % of 𝐶 1 YZY ) on the (normalized) exciton density decay. Fixed parameters are: 𝑅 1 * 1 = 5 nm; 𝑅 1 * 1 * = 4 nm; 𝑅 1 * H = 2 nm; 𝜏 H = 4 ns, 𝐶 H = 0.25 % of 𝐶 1 YZY ; 𝑛 1 (0) = 7 % of 𝐶 1 YZY .

Figure 8 .

 8 Figure 8. The effect of the quencher lifetime 𝜏 H on the (normalized) exciton density decay considering different timescales for the dye lifetime 𝜏 1 , including (A) nanosecond (𝜏 1 = 4.242 ns), (B) picosecond (𝜏 ] = 4 ps), and (C) femtosecond (𝜏 1 = 400 fs) timescales. Fixed parameters are: 𝑅 1 * 1 = 5 nm; 𝑅 1 * 1 * = 4 nm; 𝑅 1 * H = 2 nm; 𝐶 H = 0.25 % of 𝐶 1 YZY ; 𝑛 1 (0) = 7 % of 𝐶 1 YZY .

This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 5

This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 7

Acknowledgements

The acknowledgements come at the end of an article after the conclusions and before the notes and references.

Author Contributions

We strongly encourage authors to include author contributions and recommend using CRediT for standardised contribution descriptions. Please refer to our general author guidelines for more information about authorship.

Conflicts of interest

In accordance with our policy on Conflicts of interest please ensure that a conflicts of interest statement is included in your manuscript here. Please note that this statement is required for all submitted manuscripts. If no conflicts exist, please state that "There are no conflicts to declare".


G Hassan
email: ghina.hassan@lirmm.fr

A Chemori
email: ahmed.chemori@lirmm.fr

M Gouttefarde
email: marc.gouttefarde@lirmm.fr

M El

C Francis
email: clovis.francis@ensam.eu

P E Hervé

D Sallé

A New Augmented RISE Feedback Controller for Pick-and-Throw Applications with PKMs

Keywords: RISE feedback control, model-based robust control, parallel kinematic manipulator, pick-and-throw, numerical simulations

Most of the current industrial applications (like food packaging, waste sorting, machining, etc.) use parallel kinematic manipulators (PKMs) owing to their high speed and accuracy. However, parallel robots are exposed to highly nonlinear dynamics, time-varying parameters and uncertainties, especially in those applications. Considering all these issues, the synthesis of advanced and robust control schemes for PKMs is considered a challenging task. A new control scheme based on the Robust Integral of the Sign of the Error (RISE) control scheme is proposed in this work. A revision of the standard RISE control law is proposed by considering, in the control loop, a compensation term computed from the dynamic model of the robot, the measured and the desired trajectories, and the tracking error. In addition, we propose to extend the resulting controller with a nonlinear feedback function to compensate for the errors resulting from using the desired trajectories instead of the measured ones in the dynamic compensation term. The proposed control contribution can compensate for PKM parameter uncertainties and high nonlinearities as well as improve the robustness of the standard RISE controller. Numerical simulations have been conducted on a parallel robot, called T3KR, in a pickand-throw task under different operating conditions to confirm the effectiveness of the proposed control scheme.

INTRODUCTION

Parallel Kinematic Manipulators (PKMs) will, in a near future, be increasingly used in industries due to their high dynamics, rigidity, high accuracy and high payload-to-weight ratio [START_REF] Natal | Dual-space control of extremely fast parallel manipulators: Payload changes and the 100g experiment[END_REF]. A parallel robot is generally defined as follows [START_REF] Merlet | Parallel robots[END_REF]: "A parallel kinematic manipulator is a structure whose moving platform is connected to a fixed base by several independent kinematic chains". Given the aforementioned advantages, PKMs have been widely applied in various applications. For example, they have been used as a flight simulator [START_REF] Stewart | A platform with six degrees of freedom[END_REF], in the packaging industry [START_REF] Clavel | Device for the movement and positioning of an element in space[END_REF], in medical applications [START_REF] Dalvand | Remote centre-ofmotion control algorithms of 6-rrcrr parallel robot assisted surgery system (pramiss)[END_REF], in machining tasks [START_REF] Escorcia-Hernández | A new solution for machining with ra-pkms: Modelling, control and experiments[END_REF], in highspeed handling (Shang and Cong, 2009) and as haptic devices [START_REF] Grange | Overview of the delta haptic device[END_REF]. In the last few years, they have been used in selective waste sorting, alongside the conventional sorting machines, to increase the purity of streams by removing the unwanted material not selected by the machines. One example is the ABB's Delta robot used as the basis of a sorting robot adopting the pick-and-place (P&P) technique BHS (2018). Recently, the pick-and-throw (P&T) approach has been applied in waste industry aiming to increase the productivity of a PKM [START_REF] Raptopoulos | Robotic pick-and-toss facilitates urban waste sorting[END_REF]. Indeed, parallel manipulators are well known for their unfavorable nonlinearities, abundant modeling uncertainties, external disturbances, and parameter variations (i.e., payload), especially in high-speed industrial applications. In order to achieve better output trajectory tracking despite all of the above challenges, a controller must be robust enough to counteract external disturbances as well as rich in dynamic model knowledge to compensate for system uncertainties and nonlinearities. RISE feedback control, developed by Xian [START_REF] Xian | A continuous asymptotic tracking control strategy for uncertain nonlinear systems[END_REF], is a good candidate for the control of highly nonlinear uncertain dynamical systems. It is a continuous nonlinear control strategy based on the integral of the sign function in terms of the tracking error ensuring disturbances rejection. It is a non-model-based controller that can guarantee a semi-global asymptotic tracking under limited assumptions on the system and the external disturbances. RISE control strategies known by their robustness and disturbance rejection ability were successfully implemented in various real-time applications [START_REF] Feemster | Jitter reduction in a directed energy application using rise[END_REF][START_REF] Yao | Rise-based precision motion control of dc motors with continuous friction compensation[END_REF]. Moreover, the high efficiency of RISE control schemes has been proved experimentally on different robotic applications such as Autonomous Underwater Vehicle (AUV) [START_REF] Fischer | Nonlinear control of an autonomous underwater vehicle: A rise-based approach[END_REF], exoskeleton devices [START_REF] Sherwani | Rise-based adaptive control for eicosi exoskeleton to assist knee joint mobility[END_REF], rigid parallel manipulator [START_REF] Escorcia-Hernández | A new solution for machining with ra-pkms: Modelling, control and experiments[END_REF] and cable-driven parallel robots [START_REF] Hassan | Rise feedback control of cable-driven parallel robots: Design and real-time experiments[END_REF]. The RISE-based controllers show in the many researches work a high performance compared to standard RISE control law. For instance, in [START_REF] Bennehar | A new rise-based adaptive control of pkms: design, stability analysis and experiments[END_REF], a RISE-based adaptive control has been proposed as a solution for the control problem of PKMs. It consists in adding model-based adaptive feedforward term to the control loop in order to compensate for parameter variations. Moreover, a Time-Varying RISE-based Control consists in replacing the static feedback gains by nonlinear ones has been proposed and validated experimentally on a parallel manipulator [START_REF] Saied | A new time-varying feedback rise control for second-order nonlinear mimo systems: theory and experiments[END_REF]. In (Escorcia-Hernández et al., 2020a), a RISE feedback controller has been extended by an adaptive feedforward compensation term based on B-Spline Neural Networks (BSNNs) to improve tracking performance of PKMs. RISE and RISE-based control strategies show, in different applications, significant performance and robustness against disturbances and uncertainties. Therefore, the improvement of this controller attracts the attention of many researchers. The RISE feedback control law is a non-model based controller that depends only on system states, which can lead to low performance in the presence of large uncertainties and hard nonlinearities. Enriching the control loop of original RISE controller by a compensation term based on the dynamic model and the system errors has the potential to improve the performance by accommodating more uncertainties and variations in the dynamic parameters. In most industrial applications, non-model based control strategies are applied due to their simplicity and ease of implementation. However, as mentioned above, PKMs are often subject to dynamic nonlinearities, uncertainties, parameter variations, external disturbances, etc. Accordingly, non-model-based controllers may lead to poor performance and even instability when operating at critical conditions (e.g. high speed applications, payload changes). Many research works show over the last decades that enriching a controller with knowledge on the manipulator dynamics can compensate for nonlinearities and parameter uncertainties, especially for high-order nonlinear systems [START_REF] Kelly | Control of robot manipulators in joint space[END_REF][START_REF] Slotine | Applied nonlinear control[END_REF][START_REF] Ren | Experimental comparison of control approaches on trajectory tracking control of a 3-dof parallel robot[END_REF]. Therefore, different model-based controllers have proposed and implemented on PKMs, such as PD with desired gravity compensation [START_REF] Kelly | Pd control with desired gravity compensation of robotic manipulators: a review[END_REF], computed torque (CT) control which uses full knowledge of the nonlinear system dynamics [START_REF] Codourey | Dynamic modeling of parallel robots for computed-torque control implementation[END_REF] and augmented PD (APD) where the dynamic term depends on both, the desired and measured states (Shang et al., 2009). Furthermore, an extended modelbased adaptive feedforward control L 1 has been developed to control mechanical manipulators [START_REF] Bennehar | Extended model-based feedforward compensation in 1 adaptive control for mechanical manipulators: Design and experiments[END_REF]. A model predictive control (MPC) is designed and implemented on a redundantly actuated parallel robot to improve tracking capability and robustness under disturbances and parameter variations [START_REF] Wen | The study of model predictive control algorithm based on the force/position control scheme of the 5-dof redundant actuation parallel robot[END_REF]. In this paper, a revision of the original RISE feedback law is performed by augmenting its control loop with a nonlinear dynamic compensation term. This term is computed based on a combination of the dynamic parameters of the controlled system, which is a PKM in our study, its state errors, and the measured and desired trajectories. In addition, the resulting controller is extended by an auxiliary nonlinear feedback term to account for the errors that occur by replacing the actual trajectories (specifically the actual velocities and accelerations), in the dynamic compensation term, with the desired trajectories. Numerical simulations are carried out in P&T scenarios with a PKM under different operating conditions (i.e., payload changes, different speeds) to demonstrate the performance improvement and enhancement provided by the proposed control scheme. The rest of this paper is organized as follows: In Section 2, the description and modeling of T3KR PKM are presented. The background on the RISE control scheme is detailed in Section 3. Section 4 is dedicated to the proposed contribution on RISE control. The obtained numerical simulation results are presented and discussed in Section 5. Section 6 concludes this paper and presents future work.

T3KR PKM: DESCRIPTION AND MODELLING

This section describes the design of the mechanical structure of T3KR robot with its kinematic and dynamic model.

Description and kinematics of T3KR PKM

T3KR is a "Pick-and-Place" PKM, developed in collaboration between Tecnalia, LIRMM and SATT AxLR. This robot is designed to be as economical in footprint as possible (i.e. all its geometrical parameters are optimized). It consists of four driven kinematic chains connected to a common mobile platform. Each kinematic chain is a series arrangement of a revolute actuator, a reararm and a forearm (composed of two parallel rods forming a parallelogram). The forearms are connected, on their first side, to the reararms by passive spherical joints and linked on their other side to the mobile platform by means of the same aforementioned joints (as shown in Fig. 1). This design allows the end-effector to have five Degrees-of-Freedom (DOFs): three translations denoted by x, y, z and the rotation ψ pf the mobile platform around the z-axis are provided by the four actuated arms. In addition, the rotation φ of the end-effector around the z-axis is provided by an actuator fixed on the platform. It should be noted that the rotational degree of freedom ψ is kinematically redundant. Indeeed, the platform is designed such that its rotation is a parallelogram mechanism movement, where the tool control point (TCP) is on the neutral axis of this mechanism. Therefore, a rotation ψ of the mobile platform does not induce a movement of the TCP, this is why ψ is permanently maintained at zero. In this work, we are concerned only with the control of the four main actuators of T3KR to perform the pick-and-throw task. Thus, we consider the 4-dimensional coordinate vector X = [x, y, z, ψ] T as a representation of the pose of the robot's end-effector, and the 4-dimensional coordinate vector q = [q 1 , q 2 , q 3 , q 4 ] T as the actuated joint positions. The differential kinematic relationship between the Cartesian and joint velocities can be written as Ẋ = J q, where Ẋ and q are the Cartesian and joint velocities, respectively and J is the Jacobian matrix.

Dynamic modelling of T3KR PKM

Model-based controllers require the inclusion of an accurate dynamic model in the control loop. In order to exploit the dynamic modeling of T3KR robot, the following assumptions often considered in the computation of PKM dynamic model as a good compromise between precision and complexity. [START_REF] Codourey | Dynamic modeling of parallel robots for computed-torque control implementation[END_REF]: Assumption 1: The masses of the forearms are smaller than the other parts of the robot, hence their inertia is neglected. Assumption 2: The mass of each forearm is divided into two point masses located at both extremities of the forearms. Based on the virtual work principle explained in [START_REF] Bennehar | A new rise-based adaptive control of pkms: design, stability analysis and experiments[END_REF][START_REF] Escorcia-Hernández | A new solution for machining with ra-pkms: Modelling, control and experiments[END_REF], the dynamics of the T3KR robot can be reduced to the analysis of the dynamics of two main bodies: The moving platform and the main action mechanism consisting of actuators in conjunction with their corresponding reararms and forearms. Regarding the moving platform's dynamics, one can define two kinds of forces acting on it induced by gravity and Cartesian accelerations:

G t p = -M t p G, F t p = M t p Ẍ (1)
where M t p = diag{m t p , m t p , m t p , i t p } with m t p = m n + 4 m f 2 is the total mass of the moving platform with the half-masses of the forearms. i t p is the total inertia of the mobile platform. G = [0, 0, g, 0] T is the gravity vector, being g = 9.81m/s 2 the gravity acceleration, and Ẍ ∈ R 4 gives the Cartesian acceleration vector. The contributions of G t p and F t p to actuator torques can be computed by means of the Jacobian matrix J ∈ R 4×4 as follows:

Γ G t p = J T G t p , Γ F t p = J T F t p
(2) Regarding the dynamics of the actuators with their corresponding reararms and forearms, three contributing torques can be distinguished: (i) the actuators input torque Γ, (ii) the torque contribution of the gravitational forces acting on the reararms Γ G arm , and (iii) the torque contribution of the inertial force acting on the reararms Γ arm . The last two ones are denoted as follows:

Γ G arm = -gM r Cos(q), Γ arm = I arm q (3) where M r = diag{m req , m req , m req , m req }, being m req = m r L rG + L m f 2
, m r is the mass of each reararm, l rG is the distance from the axis of rotation of each reararm to its center of gravity, and L is the complete length of each reararm. Cos(q) = [cos(q 1 ), cos(q 2 ), cos(q 3 ), cos(q 4 )] where q i , for i = 1...4 are the measured joint positions. We have I arm = I act + I r + L 2 m f 2 , being given that I act and I r are the actuators inertia and the reararm inertia, respectively. The term L 2 m f 2 corresponds to the inertial contribution of the forearms using the second assumption, and q represents the accelerations in joint space. Following [START_REF] Codourey | Dynamic modeling of parallel robots for computed-torque control implementation[END_REF], stating that the sum of all noninertial forces should be equal to the sum of all inertial forces, and after rearranging the terms, the inverse dynamic model of T3KR robot can be expressed in terms of the joint coordinates q as follows: M(q) q +C(q, q) q + G(q) = Γ(t) (4) where M(q) = I arm + J T M t p J is the total mass and inertia matrix of the robot, C(q, q) q = J T M t p J denotes the Coriolis and centrifugal forces matrix, G(q) = -Γ G arm -Γ G t p represents the gravitational forces vector, Γ(t) is the control input vector. The main geometric and dynamic parameters of T3KR parallel robot are summarized in Table 1.

BACKGROUND ON RISE FEEDBACK CONTROL LAW

The dynamics of a n-DOFs kinematic manipulator can be described in joint space as follows [START_REF] Cheng | Dynamics and control of redundantly actuated parallel manipulators[END_REF]:

M(q) q +C(q, q) q + G(q) + F(q, q) + D(t) = Γ(t)

(5) where M, C and G are defined above in (4). F(q, q) ∈ R n represents the friction effects vector and D(t) ∈ R n is a general nonlinear disturbances vector (i.e. external disturbances, interaction with the environment, etc). RISE feedback law is a robust nonlinear strategy that exhibits semi-global asymptotic tracking. This technique includes, in addition to the proportional integral part, a unique integral sign function that constitutes the robustness term of RISE controller. In contrast to many existing robust controllers in the literature, RISE can generate continuous control signals, thus avoiding chattering effects and improving tracking performance. It can compensate for a large class of general uncertainties and external disturbances based on limited assumptions about the disturbances and the controlled system defined as follows [START_REF] Xian | A continuous asymptotic tracking control strategy for uncertain nonlinear systems[END_REF][START_REF] Bennehar | A new rise-based adaptive control of pkms: design, stability analysis and experiments[END_REF]: Assumption 3: The inertia matrix M(q) is a symmetric positive-definite matrix satisfying the following inequality bounds for all γ ∈ R n : m∥γ∥ 2 ≤ γ T M(q)γ ≤ m(q)∥γ∥ 2 (6) where m is a known positive constant, m:R ≥0 → R ≥0 is a nondecreasing function and ∥.∥ denotes the Euclidean norm of a vector. Assumption 4: C(q, q) and G(q) are bounded if q(t), q(t) are bounded and measurable. Moreover, the first two partial derivatives of M(q), C(q, q), G(q) and F(q, q) with respect to q(t) and q(t) exist and are bounded. Assumption 5: The disturbance vector and its first two time derivatives are bounded (i.e. D(t), Ḋ(t), D(t) ∈ L ∞ ) Assumption 6: The desired trajectory q d (t) ∈ R n is differentiable till the 4 th order and its derivatives are bounded. Let us define the desired joint positions, velocities and accelerations as q d (t), qd (t), qd (t), respectively. In order to develop the closed-loop error system, the combined velocity-position tracking error e 2 ∈ R n and the auxiliary filtered tracking error r ∈ R n are denoted as follows:

e 2 = ė + Λ 1 e. ( 7
)
where e = q dq is the output tracking error in joint space. q d ∈ R n is the desired joint positions and q is the measured one by means of the actuators encoders. Λ 1 ∈ R n×n is a diagonal, positive-definite gain matrix.

Based on the closed-loop stability conditions detailed in [START_REF] Xian | A continuous asymptotic tracking control strategy for uncertain nonlinear systems[END_REF], RISE feedback control law is given as follows:

Γ RISE = (K s + I)e 2 (t) -(K s + I)e 2 (t 0 ) + t t 0 [(K s + I)Λ 2 e 2 (σ ) + β sgn(e 2 (σ ))]dσ . ( 8
)
where K s , Λ 2 , and β ∈ R n×n are positive-definite, diagonal gain matrices, I ∈ R n is identical matrix, t 0 is the initial time and sgn is the vector of the sign functions of the combined tracking error. It is worth noting that the second term of the R.H.S of (8) (i.e., (K s + I)e 2 (t 0 )) is introduced to guarantee a zero control input at time t = t 0 (i.e., Γ(t 0 ) = 0). For more details about the stability analysis of RISE feedback control, the reader can refer to [START_REF] Xian | A continuous asymptotic tracking control strategy for uncertain nonlinear systems[END_REF].

PROPOSED CONTRIBUTION: EXTENDED RISE CONTROL PLUS COMPENSATION

General Overview on PD Control plus Compensation

Several model-based controllers have already been reported in the introduction, highlighting their improved overall tracking performance compared to non-model-based controllers. The PD plus compensation control law, called PD + , is one of the model-based controllers that is popular in academia. It is a non-adaptive version of the first adaptive controller proposed in 1987 by Slotine and Li, and it is referred to by the names of its creators: "Slotine and Li controller" [START_REF] Slotine | On the adaptive control of robot manipulators[END_REF]; [START_REF] Slotine | Adaptive manipulator control: A case study[END_REF]. From a structural point of view, this controller is composed of a PD feedback term plus a compensation term based on the full knowledge of the dynamic model combined with the system state errors. Let us first consider the classical PD control law as follows:

Γ PD = K p e(t) + K d ė(t) (9
) where K p and K d ∈ R + are constant feedback gains that are tuned to ensure the stability of the system, and e(t) = q dq is the tracking error. Enriching this controller explicitly with the dynamic model of the nonlinear controlled system allows to compensate for uncertainties and parameter variations, and thereby, enhancing the trajectory tracking accuracy. The following equation describes the control law of PD plus compensation:

Γ PD + = M(q)( qd + α ė(t)) +C(q, q)( qd +αe(t)) + G(q) + K p e(t) + K d ė(t) (10)
where α ∈ R + is defined as:

α = K -1 d K p .

Proposed Extended RISE Control plus Compensation

RISE is a non-model based controller consisting of two main parts: A linear state feedback term similar to a PI controller, depending on the combined tracking error, and a nonlinear robustness term based on the integral of the sign of the combined error. It does not take advantage of knowledge of the manipulator's dynamic model in its control loop. Thus, the dynamic parameter uncertainties and system nonlinearities are not well compensated for by the standard RISE feedback law, especially under critical operating conditions. These issues may lead to high-gain or high-frequency feedback and poor performance in the case of high nonlinearities or in presence of large disturbances. Therefore, to improve the performance of the original RISE controller, we propose to enrich its control structure by adding a full dynamic compensation term. As it can be seen from the equation of PD control plus compensation (10), the compensation term consists of the system dynamic parameters computed online using the measured trajectories and then multiplied by the sum of the desired trajectories with the state errors. To overcome the error resulting from using the desired trajectory signals, RISE control plus compensation is revised by adding a nonlinear auxiliary term depending on the tracking error e 1 and the combined error e 2 . This auxiliary term has been added to the Desired Compensation Adaptive law (DCAL) proposed in [START_REF] Sadegh | Stability and robustness analysis of a class of adaptive controllers for robotic manipulators[END_REF] to compensate for the additional error emerging from using desired trajectories instead of the measured ones in the adaptive feedforward term.

The proposed extended RISE control plus compensation, called ERISE + , is expressed as follows:

Γ ERISE + = M(q)( qd + Λ 1 ė(t)) +C(q, q)( qd + Λ 1 e(t)) + G(q) + (K s + I)e 2 (t) -(K s (t 0 ) + I)e 2 (t 0 ) + t t 0 δ ∥e(σ )∥ 2 e 2 (σ )dσ + t t 0 [(K s0 + I)Λ 2 e 2 (σ ) + β sgn(e 2 (σ ))]dσ
(11) Based on the stability analysis of the proposed controller, the nonlinear auxiliary term inspired by the DCAL scheme was added with an integral term. Due to space limitation, the stability analysis is not provided in this paper.

NUMERICAL SIMULATION RESULTS

In this section, the obtained simulation results of the PD control with computed, the original RISE, proposed ERISE + .

PD control with computed feedforward

PD control with computed feedforward consists of using the full inverse dynamic model to compensate the effect of nonlinearity but within an offline-computation mode. In fact, the inverse dynamic model is evaluated with the desired trajectories instead of the measured and estimated ones. The joint space control law can be formulated as follows [START_REF] Bennehar | Some contributions to nonlinear adaptive control of PKMs: from design to real-time experiments[END_REF]:

Γ FFPD = M(q d ) qd +C(q d , qd ) qd + G(q d ) + K p e + K d ė (12)
where K p , K d ∈ R n×n are positive definite diagonal feedback gain matrices. e(t), ė(t) ∈ R n are the joint position and velocity tracking errors, respectively.

Pick-and-Throw reference trajectory generation

The reference trajectories, illustrated in Fig. 2, are generated in Cartesian space using a third-order S-curve polynomial motion profile. These generated trajectories correspond to the scenario where the robot has to successively throw three objects of different masses to a target position, P f , located outside the robot's workspace. First, the robot has to move from the central position P 0 to the first pick position, P 1 , to grasp the detected object. Then, and according to the pick and target positions of the corresponding object, a release position P r 1 is calculated. After finding the appropriate release point, the robot accelerates to this point (i.e., P r 1 ), and then throws the object to the desired target P f . After throwing the object, the robot decelerates to pick the next object, while the released object follows its free flight ballistic trajectory from P r 1 to P f . The same cyclic movement is repeated for the second and third objects, located respectively at P 2 and P 3 . After releasing the last object, the robot returns to its initial position P 0 . The whole P&T trajectory is generated taking into account the maximum dynamic performance of the robot (i.e., maximum jerk, acceleration and velocity), the limits of its workspace and the desired target position. Referring to Fig. 2, the red lines represent the portions of the trajectory where the robot carries the object, while the green lines correspond to the portions after the release point where the robot is moving without payload. For more details on the generation of the robot's throw motion with minimum time, the reader can refer to [START_REF] Hassan | Time-optimal pick-and-throw s-curve trajectories for fast parallel robots[END_REF].

Performance evaluation criteria

To quantify the enhancement provided by the proposed ERISE + controller, we consider the Root Mean Square Error (RMSE) criterion for Cartesian positions RMSE x and joint positions RMSE q . These evaluation criteria given by the following equations allow to evaluate the effectiveness of each new control algorithm:

RMSE x = ( 1 N N ∑ i=1 (e 2 x (i) + e 2 y (i) + e 2 z (i)) (13)
RMSE q = ( 1 N N ∑ i=1 (e 2 q 1 (i) + e 2 q 2 (i) + e 2 q 3 (i) + e 2 q 4 (i)) (14)
where e x , e y and e z denote the Cartesian position tracking errors along x, y and z axes, respectively. e q 1 , e q 2 , e q 3 and e q 4 being the joint position tracking errors, and N is the number of samples.

Control gains tuning procedure

Using the trial and error method, the feedback gains of all the implemented controllers are adjusted. Manually and continuously, we edit the different sets of control gains until the desired performance is obtained. The control gains of the proposed ERISE + controller are tuned in the same way of those of the original RISE controller using the following procedure: (i) Set Λ 2 = 0, β = 0 and δ = 0, (ii) Set Λ 1 and K s as if it was a PD control plus compensation, where Λ 1 (K s + 1) is the proportional gain and (K s + 1) is the derivative gain until satisfactory tracking of the reference position and velocity signals is achieved. Then, (iii) start increasing Λ 2 by changing again Λ 1 and K s either increasing or decreasing until the best possible performance is reached. (iv) Increase β gradually to avoid chattering Table 2. Summary of the feedback control gains.

RISE FFPD Proposed ERISE + Λ 1 = 674 β = 2.5 K p = 2863 Λ 1 = 150 β = 2.5 Λ 2 = 7.2 K d = 22 Λ 2 = 1.5 K s = 21 K s = 21
effects, to improve the robustness of the controller towards disturbances. Finally, (v) increase δ to improve the overall performance while keeping control input torques below saturation. The obtained control gains for the proposed ERISE + , the standard RISE and PD with feedforward controllers are summarized in Table 2.

Obtained simulations results

To demonstrate the performance of the proposed ERISE + controller, a comparative study has been performed with the original RISE and the PD plus feedforward controller through numerical simulations on T3KR robot in a P&T task. This comparison has been established in Matlab/Simulink environment with sampling time of 0.4 ms, using the reference P&T trajectory depicted in Fig. 2. Two main scenarios have been implemented on this validation: 1) scenario 1: Robustness towards payload changes, 2) scenario 2: Robustness towards speed variations. For more realistic simulations, white noise has been added to the output joint positions, coulomb and viscous friction has been added to the dynamic model of the robot as well as a 25% of uncertainty on the inertia value I arm has been considered.

Scenario 1 -Robustness towards payload changes: In this scenario, the maximum operating acceleration was set to 4.2 G. It is the minimum sufficient value required by T3KR robot to throw an object outside of its workspace. The three objects, used for this demonstration, have different masses, which allows to evaluate the robustness of the proposed controller towards variations in payload. The first object has a mass of 50 g, the second one has a mass of 100 g (i.e. ∆ mass = +100 % w.r.t the first object), while the third one has a mass of 150 g (i.e. ∆ mass = +200 % w.r.t the first object).

The Cartesian tracking errors of the three tested controllers are plotted in Fig. 3. The obtained results clearly show the superiority of the proposed controller over the other two ones along all axes. The RMSE performance indices are evaluated for all controllers, in both Cartesian and joint spaces, and the obtained results are reported in Table 3. According to these indices, the proposed ERISE + outperforms the standard RISE by 82.7% and 72.6% for the Cartesian and joint spaces, respectively. Compared to the FFPD control, the proposed ERISE + improves the tracking performance by up to 50.3% and 45.9% in the Cartesian and joint spaces, respectively. The evolution of the generated control input torques for all controllers is displayed in Fig. 4. The control signals show, for all controllers, a good and smooth behavior within the admissible limits of the actuators of the robot (the maximum torque of T3KR actuators is 28.9 N.m). In addition, a slight reduction in energy consumption is notified for the proposed ERISE + controller compared to the two other controllers. This scenario confirms the effectiveness of the proposed controller compared to the standard RISE and FFPD controllers. The proposed ERISE + control scheme is more robust towards variations in payload, thus, it is more suitable for arbitrary P&T applications such as waste sorting.

Scenario 2 -Robustness towards speed variations: The T3KR robot is intended to be used for high-speed P&T sorting applications. Accordingly, it is useful to evaluate the tracking performance of the proposed ERISE + strategy at highacceleration conditions. The operating acceleration is increased up to 9 G. In this scenario, the robot performs the same P&T trajectory with the same manipulated objects as in the previous scenario. In Fig. 5, one can obviously see the significant improvements obtained by the proposed ERISE + control scheme along all translational axes. These improvements are quantified by exploiting the RMSE evaluation criteria in Cartesian and joint spaces. The obtained results, summarized in Table 3, show improvements of 81.4% in the Cartesian space and 71.1% in the joint space compared to the standard RISE controller. In comparison to the FFPD control law, the tracking performance is improved by up to 31.3% and 28.7% in the Cartesian and joint spaces, respectively. The evolution of the control input torques, generated by the three controllers, are depicted in Fig. 6. It can be seen that all the control signals are continuous and evolve within the admissible range of the actuators' capabilities. Moreover, as shown in Fig. 6, the proposed controller slightly reduces the power consumption as it generates less input torques, compared to the standard RISE and FFPD controllers. The overall performance improvement, obtained by the proposed ERISE + scheme, can be explained by the good compensation of the system nonlinearities provided by the contribution of the designed compensation dynamic term.

CONCLUSION AND FUTURE WORK

The main contribution of this work is a new category of RISE control strategy applied to PKMs. RISE is a non-model based

  Fig. 1. A CAD View of T3KR PKM with its main components.

  Fig. 2. 3D-view of the P&T reference trajectories of T3KR robot with the ballistic motions of the thrown objects.

Fig. 3 .

 3 Fig. 3. Scenario 1: Evolution of the Cartesian tracking errors versus time.

Fig. 4 .

 4 Fig. 4. Scenario 1: Evolution of the control input torques versus time.

Fig. 5 .

 5 Fig. 5. Scenario 2: Evolution of the Cartesian tracking errors versus time.

Fig. 6 .

 6 Fig. 6. Scenario 2: Evolution of the control input torques versus time.controller that provides semi-global asymptotic tracking of the reference trajectory under certain limitations of the system dynamics and external disturbances. We proposed to add a nonlinear model-based compensation term inspired by PD control with compensation. The compensation term consists of dynamic parameters calculated online and multiplied by the sum of the desired trajectories with state errors. To overcome the error resulting from the use of desired joint accelerations and velocities, the resulting controller is extended by a nonlinear auxiliary function based on the tracking error. The standard RISE control, PD plus feedforward and the proposed extended RISE control with compensation have been implemented through numerical simulations on the T3KR parallel robot. The obtained results clearly show the superiority of the proposed controller over the other two controllers in terms of tracking accuracy and robustness w.r.t. payload and velocity changes. This work can be extended with the stability analysis of the proposed controller as well as its validation in real-time experiments. In addition, the dynamic parameters of the manipulator may vary over time or be unknown, thus, a real-time estimation of these modeled parameters may be considered to further improve its tracking performance.

Table 1 .

 1 Summary of the main geometric and dynamic parameters of T3KR PKM.

	Parameter	Value	Parameter	Value
	L (Reararm length)	400 mm	m n (Nacelle mass)	5.68 Kg
	l (Forearm length)	900 mm	I act (Actuator inertia) 0.000969 Kg.m 2
	m r (Reararm mass)	3.28 Kg	I r (Reararm inertia)	0.173723 Kg.m 2
	m f (Forearm mass)	0.8 Kg

Table 3 .

 3 Summary of the obtained tracking performances.

	Scenario	Control	RMSE x [mm]	RMSE J [deg]
		RISE	0.0736	0.0073
	Scenario 1	FFPD Proposed ERISE +	0.0256 0.0127	0.0037 0.0020
		IMP w.r.t RISE	82.7 %	72.6 %
		IMP w.r.t FFPD	50.3 %	45.9 %
		RISE	0.1393	0.0180
	Scenario 2	FFPD Proposed ERISE +	0.0377 0.0259	0.0073 0.0052
		IMP w.r.t RISE	81.4 %	71.1 %
		IMP w.r.t FFPD	31.3 %	28.7 %

ACKNOWLEDGEMENTS

This work is supported by Tecnalia "Pick-and-Throw" research project. The corresponding author acknowledged the Lebanese University and the social foundation AZM and SAADE for financial support.