Ryan N Doan

Elaine T Lim

Silvia De Rubeis

Catalina Betancur

David J Cutler

Andreas G Chiocchetti

Lynne M Overman

Aubrie Soucy

Susanne Goetze

Christine M Freitag

Sequencing Autism

Consortium

Mark J Daly

Christopher A Walsh

Joseph D Buxbaum

Timothy W Yu
email: timothy.yu@childrens.harvard.edu

Recessive gene disruptions in autism spectrum disorder

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

Editorial Summary:

Analysis of whole exome sequencing data from 2,343 individuals with autism spectrum disorder (ASD) compared to 5,852 unaffected individuals demonstrates an excess of biallelic, autosomal mutations for both loss-of-function and damaging missense variants.

We performed a systematic analysis of exome data from the Autism Sequencing Consortium (ASC) 15 , representing 2,343 affected and 5,852 unaffected individuals (Supplementary Table 1, Supplementary Figure 1). We cataloged a total of 696,143 autosomal loss-of-function (LOF) events (representing 28,685 unique variants in 11,745 unique genes) that introduce a stop codon or disrupt a canonical splice site (Table 1). After excluding common variants (allele frequency (AF) >1%), there were 84,645 rare LOF events (27,648 unique alleles) for an average of ~10 LOF mutations per individual. After computational phasing, we found 298 events in the ASC (after filtering to exclude common polymorphisms) consistent with complete gene knockout (homozygous or compound heterozygous LOF mutations), affecting 266 individuals (Table 1). Affected individuals were disproportionately likely to harbor a gene knockout (62% more likely; 0.05 events per affected individual versus 0.031 per unaffected individual, p=0.0003 by random permutation testing), with the bulk of the excess arising from the rarest alleles (Table 1; Supplementary Table 2; Figure 1a). To control for possible differences in population and family structure (e.g., founder effects and/or consanguinity in the Finnish and Middle Eastern cohorts), we also normalized to background burdens of biallelic synonymous variants (see Methods and Supplementary Table 3-6). Individuals with ASD continued to exhibit higher knockout rates after normalization (0.042 biallelic LOF events per affected individual versus 0.031 per unaffected individual, p=0.008 by random permutation testing) (Supplementary Figure 2a). Based on the observed ascertainment differentials (AD) between affected and unaffected individuals (0.050 vs. 0.031, or 0.042 vs. 0.031 after normalization), these burdens predict a contribution of biallelic LOF alleles to ~1-2% of ASD cases (see Methods).

We considered whether these findings might be extended to incorporate the impact of missense variation. Missense variation occurs much more frequently than truncating mutations, with a small yet significant subset likely to be damaging. We evaluated potential missense contributions of two different types: all missense/nonsynonymous variants (Nsyn), and damaging missense/nonsynonymous variants (NsynD4, defined as missense changes classified as damaging by 4 or more predictive algorithms, representing 4% of all Nsyn events). Individuals with ASD exhibited an excess of biallelic events specifically involving a LOF allele in trans with a NsynD4 allele (Figure 1b), an effect which persisted after normalization to synonymous rates (Supplementary Figure 2b). Similarly, considering biallelic events involving either rare LOF or NsynD4 alleles, individuals with ASD demonstrated an excess of biallelic variation (p<10 -4 ) with or without normalization (Supplementary Figure 2c, Figure 1c). Biallelic NsynD4 alleles were also found in excess in cases (Figure 1c) but this did not persist after normalization (p=0.07, Supplementary Figure 2c).

Alleles contributing to biallelic LOF knockouts in cases were not in excess in the heterozygous state (Supplementary Figure 3). Genes impacted by biallelic mutation also scored low with respect to haploinsufficiency [START_REF] Huang | Characterising and predicting haploinsufficiency in the human genome[END_REF] and gene constraint [START_REF] Lek | Analysis of protein-coding genetic variation in 60,706 humans[END_REF] (Figure 1d-e, Supplementary Table 7), compared to known dominant and de novo ASD and ID genes [START_REF] Betancur | Etiological heterogeneity in autism spectrum disorders: more than 100 genetic and genomic disorders and still counting[END_REF] .

These data support a recessive risk model, as opposed to the additive contributions of alleles acting co-dominantly.

ASD exhibits a striking male bias, and females with ASD exhibit a disproportionate burden of de novo and inherited SNVs and CNVs, supporting a female protective/male susceptibility effect [START_REF] Jacquemont | A higher mutational burden in females supports a "female protective model" in neurodevelopmental disorders[END_REF][START_REF] Robinson | Examining and interpreting the female protective effect against autistic behavior[END_REF] . We asked whether this also applies to rare recessive mutations. Unaffected males and females did not differ in burdens of LOF knockout (Figure 2a, Supplementary Figure 4a). Females with ASD, however, exhibited a larger excess of biallelic LOF mutations (Figure 2a, Supplementary Figure 4a). Similar patterns were observed when the analysis was extended to incorporate damaging missense variants, either alone (NsynD4) (Figure 2b, Supplementary Figure 4b) or grouped with LOF mutation (LOF or NsynD4) (Figure 2c, Supplementary Figure 4c). As observed before, LOF alleles were not in excess in the heterozygous state (Supplementary Figure 3). From these patterns we conclude that the female protective/male susceptibility effect extends to recessive mutations as well.

The ascertainment differentials observed led us to estimate the total contribution of biallelic mutation to ASD to be 3-5% (2% from LOF mutations, 1-3% from missense mutations).

Furthermore, based on the ascertainment differentials, identifiable biallelic LOF or damaging missense alleles appeared to contribute to ~10% of female ASD cases in the ASC (Supplementary Figure 4b), suggesting that the category of biallelic mutation could be one of the largest genetic contributors to ASD in females.

To address potential biological implications of these results, we examined the genes in which biallelic mutations were found. We first focused on the category of strict biallelic LOF mutations, estimating from the observed ascertainment differentials that ~25-40% of those found in ASD patients (0.019 out of 0.05; or 0.011 out of 0.042 after normalization to synonymous burdens) contribute to their condition. A total of 109 unique genes were knocked out (biallelic LOF) in affected individuals in our dataset. We reasoned that the most relevant of these genes would be those that are selected against biallelic inactivation in the general population. We therefore developed a filter to remove any gene that is also knocked out in one or more of 60,706 individuals sequenced as a part of ExAC (see Methods) [START_REF] Lek | Analysis of protein-coding genetic variation in 60,706 humans[END_REF] . To test its effectiveness, we examined its impact on observed patterns and burdens of gene knockouts in ASD. Overall, 131 of the 227 (58%) gene knockouts observed in the ASC were also observed in ExAC. Application of this filter reduced the overall burdens of observed biallelic knockout by half but left the ascertainment differentials intact (Figure 3a, Supplementary Figure 5). This significantly increased the observed proportional excess of biallelic LOF mutation: after filtering, on average, individuals with ASD were ~2.6 times as likely to have biallelic LOF mutations as individuals without ASD (p=0.002) (3.9X in females, p=0.001, and 2.6X in males, p=0.01) (Figure 3a).

We applied a similar logic to biallelic missense mutations as well. We reasoned that missense alleles encountered in the homozygous state in the general population are not likely to contribute to recessive condition. ~50% of rare, biallelic damaging missense events in the ASC involved an allele encountered in the homozygous state in ExAC (Supplementary Figure 6a, Supplementary Table 8). Filtering these events allowed us to discern for the first time a significant excess of rare biallelic damaging missense variation (i.e. NsynD4) in cases compared to controls (Supplementary Figure 6a), with an ascertainment differential of ~4%. Biallelic Nsyn variants were also in slight excess (p=0.0037), which arose almost entirely from affected females (Supplementary Figure 6b).

These results indicate that the observed excess in ASD is driven by genes and alleles that are biallelically constrained in the control population, strongly supporting the biological relevance of these findings.

De novo gene disrupting mutations contribute strongly to ASD [START_REF] Rubeis | Synaptic, transcriptional and chromatin genes disrupted in autism[END_REF] . To examine possible interactions between de novo and biallelic gene disrupting events, we compared the burden of biallelic LOF events in cases with and cases without LOF de novo mutations (Supplementary Figure 10) [START_REF] Rubeis | Synaptic, transcriptional and chromatin genes disrupted in autism[END_REF] . An excess of rare biallelic LOF events is strongly evident in ASD cases lacking de novo LOF variants (Supplementary Figure 10, p=0.0006). Their contribution in ASD cases with de novo LOF variants is harder to distinguish, as the observed difference in burdens of biallelic LOF events did not achieve statistical significance (Supplementary Figure 7). These may represent orthogonal risk categories, although larger sample sizes will be required to clarify this in the future.

After ExAC filtering, a total of 57 biallelically constrained genes harbored biallelic LOF mutation in at least one ASD case. Of these, we can estimate that ~60% of these (based on the observed burdens, see Methods), and ~80% of the knockouts observed in females, likely confer ASD risk. To further enrich for relevant genes, we also set aside 13 genes knocked out in one or more unaffected individuals in our ASC dataset (Figure 3b, Supplementary Table 7), and three that failed Sanger confirmation, leaving a final set of 41 genes that were specifically knocked out in individuals diagnosed with ASD, but never in 66,557 controls (60,706 from ExAC and 5,852 unaffected individuals from ASC) (Supplementary Table 9,10).

Several of these represent genes already implicated in disease. Eight (CA2, DDHD1, NSUN2, PAH, RARB, ROGDI, SLC1A1, USH2A) were associated with recessive genetic conditions involving neurodevelopmental delay (Table 2 and Supplemental Information for clinical details). In at least six of these cases (CA2, NSUN2, PAH, RARB, ROGDI, USH2A), available medical and phenotypic records were concordant with the clinical features predicted from the genetic finding (Supplemental Information). Phenotypic data was either unavailable or incomplete, but not discordant, for other cases. In two cases, knockouts were encountered in genes that cause well-characterized autosomal recessive conditions without clear connection to autism (RFX5 and DNAI2, associated with autosomal recessive immunodeficiency and primary ciliary dyskinesia, respectively). In two other cases, knockouts were observed in genes associated with autosomal dominant disorders via recurrent gain-of-function mutations (IFITM5 and KIF22, causing osteogenesis imperfecta type V and spondyloepimetaphyseal dysplasia with joint laxity type 2, respectively). The phenotypic consequence of recessive loss-of-function mutations for these two genes is unknown. Interestingly, KIF22 is one of 25 genes in the 593kb core interval affected by 16p11.2 deletions. [START_REF] Weiss | A recurrent genetic cause of autism: microdeletion at 16p11[END_REF] KIF22 was proposed as a one of two genes in the interval that demonstrates dosage sensitivity, with knockdown causing axon tract and behavioral defects in zebrafish, [START_REF] Blaker-Lee | Zebrafish homologs of genes within 16p11.2, a genomic region associated with brain disorders, are active during brain development, and include two deletion dosage sensor genes[END_REF] and disruption a KIF22 ortholog in Drosophila causes synapse development defects [START_REF] Park | Drosophila Homolog of Human KIF22 at the Autism-Linked 16p11.2 Loci Influences Synaptic Connectivity at Larval Neuromuscular Junctions[END_REF] .

Several of the remaining 33 genes are of strong neurobiological interest. One affected individual bore a homozygous LOF mutation (p.E48X) in SLC1A1, mutations in which have been linked to dicarboxylic aminoaciduria, [START_REF] Bailey | Loss-of-function mutations in the glutamate transporter SLC1A1 cause human dicarboxylic aminoaciduria[END_REF] an extremely rare disorder with prior associations with intellectual disability [START_REF] Teijema | Dicarboxylic aminoaciduria: an inborn error of glutamate and aspartate transport with metabolic implications, in combination with a hyperprolinemia[END_REF][START_REF] Swarna | Dicarboxylic aminoaciduria associated with mental retardation[END_REF] and obsessive compulsive disorder (OCD) [START_REF] Bailey | Loss-of-function mutations in the glutamate transporter SLC1A1 cause human dicarboxylic aminoaciduria[END_REF] . SLC1A1 encodes a brain-expressed, activity-regulated glutamate transporter [START_REF] Rothstein | Localization of neuronal and glial glutamate transporters[END_REF][START_REF] Ross | MB mRNA for the EAAC1 subtype of glutamate transporter is present in neuronal dendrites in vitro and dramatically increases in vivo after a seizure[END_REF] that modulates glutamatergic neurotransmission via a variety of mechanisms [START_REF] Nieoullon | The neuronal excitatory amino acid transporter EAAC1/EAAT3: does it represent a major actor at the brain excitatory synapse[END_REF][START_REF] Bianchi | Changes in the expression of the glutamate transporter EAAT3/EAAC1 in health and disease[END_REF][START_REF] Stafford | Glutamate spillover augments GABA synthesis and release from axodendritic synapses in rat hippocampus[END_REF][START_REF] Mathews | Neuronal glutamate uptake Contributes to GABA synthesis and inhibitory synaptic strength[END_REF][START_REF] Scimemi | Neuronal Transporters Regulate Glutamate Clearance, NMDA Receptor Activation, and Synaptic Plasticity in the Hippocampus[END_REF][START_REF] Diamond | Neuronal glutamate transporters limit activation of NMDA receptors by neurotransmitter spillover on CA1 pyramidal cells[END_REF] SLC1A1 knockout in mice results in altered locomotor behavior [START_REF] Peghini | Glutamate transporter EAAC-1-deficient mice develop dicarboxylic aminoaciduria and behavioral abnormalities but no neurodegeneration[END_REF] and learning and memory defects. [START_REF] Lee | Effects of isoflurane on learning and memory functions of wild-type and glutamate transporter type 3 knockout mice[END_REF] Two brothers with ASD bore homozygous stop-gain mutations in FEV, a transcription factor that is required for both development and function of serotonergic neurons (Supplementary Figure 8a). These neurons send projections that ramify broadly throughout the CNS (Supplementary Figure 8b), and have wide-ranging neuromodulatory effects on many physiological states and behaviors, including mood, aggression, anxiety, sleep, and movement [START_REF] Lucki | The spectrum of behaviors influenced by serotonin[END_REF] . Serotonergic dysfunction has been strongly implicated in neuropsychiatric conditions including ASD [START_REF] Muller | The serotonin system in autism spectrum disorder: From biomarker to animal models[END_REF] , schizophrenia and mood disorders [START_REF] Lucki | The spectrum of behaviors influenced by serotonin[END_REF] . In mice, Fev/Pet-1 has been shown to be highly expressed in the serotonergic raphe nuclei [START_REF] Hendricks | The ETS domain factor Pet-1 is an early and precise marker of central serotonin neurons and interacts with a conserved element in serotonergic genes[END_REF] , and inactivation of Fev during development [START_REF] Hendricks | Pet-1 ETS gene plays a critical role in 5-HT neuron development and is required for normal anxiety-like and aggressive behavior[END_REF] or adulthood [START_REF] Liu | Pet-1 is required across different stages of life to regulate serotonergic function[END_REF] results in depletion of serotonergic neurons and neurobehavioral changes including anxiety and aggression. [START_REF] Hendricks | Pet-1 ETS gene plays a critical role in 5-HT neuron development and is required for normal anxiety-like and aggressive behavior[END_REF][START_REF] Liu | Pet-1 is required across different stages of life to regulate serotonergic function[END_REF] Human FEV is also expressed in the serotonergic raphe in adult brain [START_REF] Iyo | Regional distribution and cellular localization of the ETS-domain transcription factor, FEV, mRNA in the human postmortem brain[END_REF][START_REF] Maurer | The Ets transcription factor Fev is specifically expressed in the human central serotonergic neurons[END_REF] , and we found similar results via in situ hybridization analysis of human fetal brain as well (Supplementary Figure 8c-f). The older affected brother was diagnosed with autism and had a measured IQ of 69. He exhibited severe, stereotyped aggressive and self-injurious behaviors. He had a BMI of 25.6 (99%) and mildly dysmorphic facial features. EEG showed bilateral rolandic focus with generalization. The younger affected brother was diagnosed with PDD-NOS and intellectual disability (IQ not measured), had a BMI of 17.7 (91%), and no facial dysmorphisms (further details in Supplemental Information). We are not aware of any previous report of a human loss-offunction phenotype for FEV.

Gene knockouts were also identified in FCHSD2, a cytoskeletal adaptor protein that is highly expressed in human fetal brain and is the human homolog of Drosophila nervous wreck [START_REF] Cao | FCHSD1 and FCHSD2 Are Expressed in Hair Cell Stereocilia and Cuticular Plate and Regulate Actin Polymerization In Vitro[END_REF] , NCOA7, a brain-specific transcriptional coactivator for steroid and nuclear hormone receptors [START_REF] Shao | ERAP140, a conserved tissue-specific nuclear receptor coactivator[END_REF] , and MOB1B, a downstream effector in the Hippo signaling pathway, a pathway recently implicated in regulation of brain size [START_REF] Poon | The Hippo Pathway Regulates Neuroblasts and Brain Size in Drosophila melanogaster[END_REF] . Many additional genes in the candidate set have abundant brain expression (BRMS1L, ELOF1, GCN1L1, NUPR1, SLC22A6) but less clear neurobiological activity.

We also expanded our scope to incorporate biallelic missense mutation. There were 409 distinct genes that harbored biallelic, damaging missense mutations in cases but not in controls (Supplementary Figure 9, Supplementary Table 11). Only one gene was found that was independently hit in multiple families: AMT, an established cause of nonketotic hyperglycinemia (NKH). Hypomorphic mutations in AMT have previously been reported to cause ASD without NKH's characteristic metabolic abnormalities. No other gene was hit more than once, so larger sample sizes would be necessary to establish any of these genes as ASD genes on this basis alone (Supplementary Table 11, Supplementary Figures 10 and11). However, there were thirteen cases in which affected individuals bore biallelic mutations in alleles of established medical significance (in ClinVar as Pathogenic or Likely Pathogenic) (Table 3, Supplementary Table 12, Supplementary Information). All were diagnostic of genetic disorders with known neurodevelopmental consequences including classic metabolic diseases, mitochondrial depletion syndrome, and other syndromic conditions. All together, combined with the previously described eight LOF cases, biallelic mutations in known genes constituted ~1% of our cohort (21 out of 2,343 affected individuals), underscoring the importance of clinical screening for monogenic recessive conditions in this patient population.

Our study extends prior estimates of the contribution of recessive mutations to ASD [START_REF] Yu | Using Whole-Exome Sequencing to Identify Inherited Causes of Autism[END_REF][START_REF] Lim | Rare Complete Knockouts in Humans: Population Distribution and Significant Role in[END_REF] and proves that the female protective / male susceptibility effect in ASD applies to biallelic mutations as well. Expanding these efforts to larger ASD cohorts should offer a path forward to delineate novel neurobiological mechanisms in autism.

METHODS

Whole exome data analyses

Whole-exome sequencing data were gathered as part of the Autism Sequencing Consortium [START_REF] Cheng | Investigating embryonic expression patterns and evolution of AHI1 and CEP290 genes, implicated in Joubert syndrome[END_REF] . We focused our analysis on nine datasets generated using the Illumina sequencing platform, Supplementary Table 1). Alignment and variant calling were performed as previously described. Low-quality variants were filtered with the following criteria: (1) SNPs with GQ<20, (2) SNPs not passing the standard PASS filter, and (3) all indels (due to excessive false positives leading to spurious frameshift calls). Variant annotation was performed using ANNOVAR [START_REF] Wang | functional annotation of genetic variants from highthroughput sequencing data[END_REF] to add allele frequencies, functional predictions, conservation, and gene annotations including OMIM and HGMD [START_REF] Stenson | The Human Gene Mutation Database: towards a comprehensive repository of inherited mutation data for medical research, genetic diagnosis and next-generation sequencing studies[END_REF] disease associations. Allele frequencies were also calculated for each contributing cohort and across the entire collection of ASC samples analyzed. Variants were categorized by their predicted impact into synonymous, nonsynonymous (missense), or loss-of-function (altering a canonical splice site or resulting in a stopgain). Additional classification of damaging missense variants (NsynD4) was performed for variants predicted to be deleterious by at least 4/6 algorithms (SIFT 51 , PolyPhen2_HDIV 52 , PolyPhen2_HVAR, MutationTaster [START_REF] Schwarz | MutationTaster evaluates disease-causing potential of sequence alterations[END_REF] , MutationAssessor [START_REF] Reva | Predicting the functional impact of protein mutations: application to cancer genomics[END_REF] , and LRT [START_REF] Chun | Identification of deleterious mutations within three human genomes[END_REF] ).

Phasing was performed using Beagle 4 [START_REF] Browning | Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering[END_REF] on the ASC dataset. VCFs were subdivided by chromosome and cohort prior to processing with Beagle employing default parameters.

Biallelic burden calculations

Allele frequency filtration was performed as described in the text using available public databases (1000G, EVS, ExAC). For rare variant calculations at maxAF 0.1%, 0.5%, 1%, 5%, and 10%, additional cohort-specific filters (cohort-specific AF of 5% or 10%) were applied to account for allelic variants that may be prevalent in a population but are not well represented in existing control databases. Burdens of homozygous and compound heterozygous variants were then counted for each cohort and for the entire ASC, stratifying by affected/unaffected status, gender, and mutational category (LOF, Nsyn, NsynD4, Syn).

To control for potential differences in underlying family and population structure of some of the constituent cohorts, we also measured rates of synonymous biallelic variation. In two cohorts (from the Middle East and Finland, respectively), rates of biallelic synonymous variation in affected individuals were higher than rates of biallelic synonymous variation in unaffected individuals (Supplementary Table 3). Therefore, for these cohorts, we normalized our LOF, NsynD4, and Nsyn burden calculations to the background rate of biallelic synonymous variation, to isolate the impact of biallelic LOF variation independent of background homozygosity in these cohorts.

The contribution of mutations to ASD risk was estimated using the following approach. First the ascertainment differential (AD) is calculated as the difference between the burden of biallelic mutations in cases and in controls and represents the estimated fraction of mutations that contribute to ASD (for example, an AD of 0.01 suggests that 1% of cases have a biallelic mutation contributing to ASD risk. The fraction of observed biallelic mutations contributing to disease is then calculated as AD divided by the total mutation burden in cases. For example, given a burden of 0.025 biallelic LOF events in cases and 0.01 biallelic LOF events in controls, the AD is 0.015 and the percentage of mutations identified which contribute to disease risk is 0.015 divided by 0.025, or 60%.

For burden comparisons, significance was assessed by permutation testing using 10,000 random assignments of diagnostic categorization to individuals within the study. A significance threshold of p<0.05 (i.e., without correction for multiple testing) was employed since the allele frequency cutoffs and variant annotation classes (LOF, LOF+NSynD4, LOF or NSynD4, NSynD4, NSyn) analyzed are not discrete. We do note that, all ExAC-filtered biallelic burden results (i.e., excess biallelic LOF, NSynD4, and NSyn mutations, summarized in Figure 3 and Supplementary Figures 5 and6) were significant even with application of maximally conservative Bonferroni correction (e.g., p<0.01 assuming five fully independent hypotheses).

ExAC filtration

To generate a list of genes commonly inactivated by biallelic LOF mutation, we analyzed whole exome data from 60,706 individuals compiled by the Exome Aggregation Consortium (ExAC) [START_REF] Lek | Analysis of protein-coding genetic variation in 60,706 humans[END_REF] . We applied ANNOVAR to systematically re-annotate all variants from ExAC against RefSeq complete protein coding genes. Since phasing was unavailable, we focused on genes bearing putative homozygous LOF variants (any variant in stopgain, frameshift, or canonical splicing change). Homopolymer expansions / contractions in canonical splice regions were ignored, as expansions and contractions in these regions are common and have little impact on splicing. With these parameters, a total of 2986 genes were found to be inactivated by homozygous stopgain, frameshift, or canonical splice site mutation in ExAC. We also used this re-annotated dataset to generate a list of ExAC alleles encountered in the homozygous state, to identify missense and other mutations that are more likely to be clinically benign. For the final LOF candidate genes presented in Supplementary Tables 9 and10, contributory alleles were also confirmed to be absent in the homozygous state in the gnomAD dataset [START_REF] Lek | Analysis of protein-coding genetic variation in 60,706 humans[END_REF] .

Variant confirmation

Sanger sequencing was employed to confirm the accuracy of exome variant calls for a subset of LOF and NsynD4 variants with an empirical validation rate of 96% (107 of 114 alleles tested). These included all LOF alleles reported here for which DNA samples were available. If DNA samples were not available, we note the presence of allelic variants in known databases and/or whether predicted genotypes follow Mendelian segregation patterns in parental exome data (where available) as additional support.

Human tissue collection and in situ hybridization for FEV

Human embryonic and fetal tissues were obtained from the MRC-Wellcome Trust Human Developmental Biology Resource (www.hdbr.org) [START_REF] Gerrelli | Enabling research with human embryonic and fetal tissue resources[END_REF] following appropriate consent from the women involved and adhering to the relevant HTA guidelines with approval from Research Ethics Committee NRES Committee North East -Newcastle & North Tyneside 1. Table 3.

Biallelic mutations identified in cases and controls previously reported as pathogenic or likely pathogenic in ClinVar.

Gene

Figure 1 .

 1 Figure 1. An excess of rare, damaging, biallelic mutation in ASD. (a) Rates of biallelic gene knockout (strict LOF) in the ASC, stratified by diagnosis and allele frequency. Rates of biallelic variation, considering (b) LOF variants paired with a damaging missense variant (NsynD4, predicted to be deleterious by at least 4 algorithms), (c) LOF or NsynD4 variants, or NsynD4 variants alone. Genes impacted by biallelic variation in cases exhibit (d) haploinsufficiency and (e) pLI score profiles consistent with recessive genes. Error bars represent the 95% confidence intervals.

Figure 2 .

 2 Figure 2. Biallelic mutations in ASD: effects of sex. (a) Rates of biallelic gene knockout (strict LOF) in the ASC, stratified by diagnosis and sex. (b, c) Rates of biallelic variation stratified by sex, considering (b) LOF or damaging missense (LOF or NsynD4) variants or (c) damaging missense variants alone (NsynD4). Error bars represent the 95% confidence intervals.

Figure 3 .

 3 Figure 3. Biallelic mutations in ASD: ExAC filtration. (a) Rates of biallelic gene knockout (strict LOF) after filtration of commonly inactivated genes in ExAC. (b) Breakdown of candidate ASD genes excluded by presence of knockouts in ExAC, or by knockouts that are common and/or in unaffected individuals of the ASC. Error bars represent the 95% confidence intervals.

  . Author manuscript; available in PMC 2019 December 17.

Table 1 .

 1 Patterns of biallelic LOF mutation in the Autism Sequencing Consortium

	Events

per individual Allele frequency # of biallelic LOF events # of unique genes involved # of individuals harboring Unaffected Affected p-value

	0.2	0.03	0.0005	0.003	<0.0001	<0.0001
	29.128	0.219	0.096	0.050	0.042	0.032
	28.997	0.194	0.068	0.031	0.022	0.016
	8195	1420	559	266	201	146
	636	414	324	237	182	138
	237,936	1,646	621	298	225	170
	All	≤10%	≤5%	≤1%	≤0.5%	≤0.1%
				Nat Genet. Author manuscript; available in PMC 2019 December 17.

Table 2 .

 2 Clinically relevant gene knockouts in the ASC

	Gene	Mutations	Disease relevance	Notes
	CA2	hom c.232+1G>A (NM_000067)	Carbonic anhydrase deficiency	Intellectual disability, osteopetrosis, renal tubular acidosis
	DDHD1	comp het c.1311-2A>T/c.2459-2A>T (NM_001160147)	Spastic paraplegia	Spastic paraplegia, mitochondrial abnormalities; only 3 cases reported
	NSUN2	hom c.C1708T:p.Q570X (NM_001193455)	Autosomal recessive intellectual disability	Syndromic intellectual disability; only 7 mutations reported
	PAH	hom c.C703T:p.Q235X and hom c. 592_613del:p.Y198fs (NM_000277) (two separate individuals)	Phenylketonuria	~20% of individuals with PKU have autism (finding previously reported 13 )
				Elevated urinary glutamate and
	SLC1A1	hom c.G142T:p.E48X (NM_004170)	Dicarboxylic aminoaciduria	aspartate, variable neuropsychiatric
				symptoms; only 3 cases reported
				Microphthalmia, pulmonary
	RARB	hom c.C78A:p.C26X (NM_000965)	PDAC syndrome	hypoplasia, diaphragmatic hernia, and cardiac defects; milder forms
				described
	ROGDI	hom c.201-1G>T (NM_024589)	Kohlschütter-Tönz syndrome	Global developmental delay, epilepsy, spasticity, amelogenesis imperfecta
	USH2A	comp het c.T12714G:p.Y4238X/ c.G6224A:p.W2075X (NM_206933)	Usher syndrome	Sensorineural hearing deficiencies, retinitis pigmentosa (finding previously reported 14 )

Nat Genet. Author manuscript; available in PMC 2019 December 17.

ACKNOWLEDGEMENTS

We thank A. Hossain and N. Hatem for their help with sample preparation, and Janet Kerwin for her help with analysis of in situ hybridization results. R.N.D. was supported by an NIH T32 fellowship from the Fundamental Neurobiology Training Grant (5 T32 NS007484-14) and the Nancy Lurie Marks Clinical and Research Fellowship Program in Autism. The ASC is supported by the National Institute of Mental Health (grants MH100233, MH100229, MH100209, MH100239, MH111661, MH111660, MH111662, and MH111658). Collection of the PAGES cohort is supported by the National Institute of Mental Health (grant MH097849). This work was supported in part through the computational resources provided by the Department of Scientific Computing at the Icahn School of Medicine at Mount Sinai, and the Research Information Technology Group at Harvard Medical School, which is partially supported by NIH grant NCRR 1S10RR028832-01. Human embryonic and fetal material was provided by the joint MRC/Wellcome Trust (grant # MR/R006237/1) Human Developmental Biology Resource (www.hdbr.org). C.A.W. is an Investigator of the Howard Hughes Medical Institute. C.A.W. and T.W.Y. were supported by the National Institute of Mental Health (NIMH) (grants RC2MH089952 and RO1MH083565). T.W.Y. was supported by NIH/NIMH R01MH113761, NICHD/NHGRI/NIH U19HD077671, NIH/NICHD U24HD0938487 and a SFARI Pilot Research Award. S.D.R. and J.D.B. are supported by the Beatrice and Samuel A. Seaver Foundation.

.v4.p3.

RNA probes were produced by amplifying human embryonic cDNA using a forward primer to FEV (NM_017521, nt. 1286-1305) tagged with a T7 RNA promoter sequence and a reverse primer (NM_017521, nt. 1862-1843) tagged with a SP6 RNA promoter sequence (Supplementary Table 13). The resulting fragment was in vitro transcribed using the DIG RNA Labeling Kit (Roche).

The in situ hybridizations were performed as previously described [START_REF] Cheng | Investigating embryonic expression patterns and evolution of AHI1 and CEP290 genes, implicated in Joubert syndrome[END_REF] , except that 60ng RNA probe in 100ul ULTRAhyb (Ambion) per slide was used for the hybridization and the posthybridization washes were modified: slides were washed for 10 minutes in 5 X SSC, twice in 2 X SSC, 50% formamide/1 X SSC, 1 X SSC and finally in 0.1X SSC. The hybridization incubation and all post-hybridization steps were performed at 63.5°C. Three human fetal brain samples were used: a Carnegie Stage (CS) 23 brain (12959), a CS23 head (12011), and a 10 post conception week hindbrain (sample 12721). Findings described were confirmed in at least two specimens. Expression patterns were analyzed relative to spatial models provided by the HuDSeN human gene expression spatial database project (http:// www.hudsen.org). [START_REF] Kerwin | The HUDSEN Atlas: a three-dimensional (3D) spatial framework for studying gene expression in the developing human brain[END_REF]

Reporting Summary

Further information on research design is available in the Life Sciences Reporting Summary linked to this article.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.

Author Manuscript

Doan et al.

Author Manuscript


Abdenbi Mohand Ousaid

Guillaume Millet

Stephane Regnier

Sinan Haliyo

Vincent Hayward
email: vincent.hayward@isir.upmc.fr

Abdenbi Mohand-Ousaid

Stéphane Régnier

Haptic interface transparency achieved through viscous coupling

Keywords: Electromagnetic drives, discretedomain pole placement, haptic interfaces

Introduction

The ideal haptic device has no mass, infinite bandwidth and can supply unlimited force or torque. In haptic device design, the actuator, or 'prime mover', is often a limiting element, followed by mechanical transmissions and structural elements which all restrict the range of signals that can be transmitted to a user [START_REF] Morrell | In pursuit of dynamic range: Using parallel coupled actuators to overcome hardware limitations[END_REF][START_REF] Hayward | Performance measures for haptic interfaces[END_REF].

Given a particular actuator, no matter how well transmissions and linkages are designed, these elements can only degrade key aspects of performance that include the dynamic range, i.e. the ratio of the largest to the smallest specifiable force, the stiffness linking the prime mover to the manipulandum, or the end-point inertia. In order to cope with these practi-cal limitations, human sensorimotor performance can be used as a factor to bound the desired performance.

One aspect of human performance that is relevant to device design is the smallest detectable force, since, ideally, this quantity should match the smallest force than can be commanded by the device. At the upper end of the range, the largest force should match the motoric output of the hand. Ideally, the dynamic range of the device should be achievable over the whole frequency range available to touch, that is, from dc to, say, 1 kHz.

While actuator saturation (short and long term) is the factor that determines the upper limit of the range, the lower limit merits discussion. When a manipulandum interacts with a hand, it is subject to several forces.1 Neglecting the effect of internal elastic forces, that is, if the device operates below its first resonant mode, these forces are: a, the actuating force, typically a Laplace force developed in the motor windings; b, forces due to induced currents which oppose motion;

c, forces due to mechanical losses: viscosity, friction; d, inertial forces to which moving parts are subjected;

e, the force applied by the hand which include at least an inertial component, a viscous component, and an elastic component, due to the movement of tissues.

According to Newton's second law, all these forces must balance. In the above list, it is well known to haptic device designers that all these factors play a significant role in the total result. To illustrate the tradeoffs involved, direct-drive, friction-free floater designs pay the price for nearly absent losses by large end-point inertias due to necessity to actuate against gravity [START_REF] Berkelman | Lorentz magnetic levitation for haptic interaction: Device design, performance, and integration with physical simulations[END_REF].

We therefore desire a transducer able to command a desired force at the interface between a manipulandum and the soft load represented by a finger that can minimize all parasitic forces. Informally, an ideally transparent device should be such that a ≈ e and that b, c, and d in the above list can be neglected. An estimate, developed further later, shows that the end-point inertia should be on the order of grams, not hundreds of grams, which is the standard in today's devices. We also desire that this force be accurately produced, not only at dc, but over a wide frequency range, up to one kilohertz.

In this article, we first develop target performance figures that result from the above considerations, and propose a dual-stage drive architecture that has the potential to meet these targets and we describe the design of its key elements. We then describe a control design performed in the sampled-data domain in order to maximize accuracy. Lastly, we tested human performance gained when using this device. To this end, we asked participants to perform a foregroundbackground signal detection task and compare the results with those obtained with a single-stage design.

Haptic Transparency

Achieving haptic transparency depends on reducing the parasitic forces under the smallest humandetectable force under all desired operating conditions. Seen from the viewpoint of the hand, we desire the device to have a mechanical impedance (the ratio of force over displacement and derivatives) that is significantly smaller than that of the fingertips since then, the deformation of the tissues would not be the result of parasitic forces, but predominantly from commanded forces. Some target values were suggested by [START_REF] Millet | The ultimate haptic device: First step[END_REF] to specify the characteristics of a haptic interface that can operate at the limits of human performance, see also [START_REF] Baud-Bovy | Hand-held object force direction identification thresholds at rest and during movement[END_REF]. These numbers, reported in Table 1, are arrived at by considering that when the hand actively interacts with the handle, parasitic forces should be at threshold, thus not interfering with the Laplace force produced by the motors.

It is worth noting that the dynamic range from the lowest to the highest specifiable force covers four orders of magnitude and that the achievable dynamic range of a high-quality dc motor is only about two orders of magnitude. It is clear that a single standard electric motor, the actuator of choice when building impedance-based haptic devices, cannot meet our objectives.

Table 1: Target performance figures with a 70 mm diameter handle [START_REF] Millet | The ultimate haptic device: First step[END_REF].

Quantity

Value Equiv. angular Value Mass a 0.5 10 -2 kg 6.1 10 -6 kg•m 2 Viscosity b 0.5 10 -1 N•s/m 6.1 10 -5 N•m•s/rad Friction c

1.0 10 -3 N 3.5 10

Related Approaches

Perhaps the oldest motivation for dual-motor designs is the desire to modulate the intrinsic dynamics of an actuated joint to imitate the antagonist action of skeletal muscles. [START_REF] Sugano | Force control of the robot finger joint equipped with mechanical compliance adjuster[END_REF]; [START_REF] English | Implementation of variable joint stiffness through antagonist actuation using rolamite springs[END_REF]; [START_REF] Stramigioli | Variable spatial springs for robot control applications[END_REF] describe efforts aimed at modulating elasticity with electric motors, and [START_REF] Boulet | Modeling and control of a robotic joint with inparallel redundant actuators[END_REF] describes the modulation of viscosity with hydraulic motors. These ideas have evolved to be included in the design of robots that can interact with people, see [START_REF] Bicchi | Safety for physical human-robot interaction[END_REF].

In the symbolism of Fig. 1, a rectangular box represents an inertia, two parallel lines represent a controllable source of force, typically a Laplace force, a coil represents a spring, a piston represents a damper, a small circle represents the load, and a mechanical connection represents a common velocity on each side. Lack of space prevents us from commenting on the various control options afforded by each design arrangement. The first case, a, is important because it models most actual haptic devices, since the preferred method of construction calls for grounding the actuators. There is an elastic element between the motor and the handle. If elasticity is built-in by design, we have the so-called "series elastic actuators" intended to be employed in closed-loop [START_REF] Pratt | Series elastic actuators[END_REF]Lauria et al., 2008). In haptics, one would rather make the transmission as stiff as possible to raise the system's natural frequency.

The next configuration, b, corresponds to the standard, variously called coarse-fine or mini-macro manipulator design [START_REF] Sharon | Enhancement of robot accuracy using endpoint feedback and a macro-micro manipulator system[END_REF][START_REF] Sharon | High bandwidth force regulation and inertia reduction using a macro/micro manipulator system[END_REF][START_REF] Salcudean | On the control of redundant coarse-fine manipulators[END_REF][START_REF] Khatib | Reduced effective inertia in macro-/mini-manipulator systems[END_REF][START_REF] Hollis | A six-degree-of-freedom magnetically levitated variable compliance fine-motion wrist: Design, modeling, and control[END_REF]. Here, a large motor 'carries' a smaller one, but both share the same load since they are on the same load path. Such a scheme has been advocated for haptic devices [START_REF] Stocco | A coarse-fine approach to force-reflecting hand controller design[END_REF][START_REF] Wall | A high bandwidth interface for haptic human computer interaction[END_REF]. With proper control, this configuration can reduce the apparent inertia of the whole system but the small motor must be able to bear the whole brunt of the load. Addressing this limitation naturally leads to an arrangement, c, described by [START_REF] Morrell | In pursuit of dynamic range: Using parallel coupled actuators to overcome hardware limitations[END_REF], and applied to making human-friendly robots and force feedback devices [START_REF] Zinn | A new actuation approach for human friendly robot design[END_REF][START_REF] Zinn | Large workspace haptic devices -a new actuation approach[END_REF].

The large motor supplies the largest portion of the force and the small motor "fills-in" during transients. With proper control, the user is exposed to the inertia of the small motor but not that of the large. A variant design is proposed in [START_REF] Conti | A hybrid actuation approach for haptic devices[END_REF], where the large motor is replaced by a passive brake.

More generally speaking, the haptic designer could use all three types of mechanical coupling: elastic, dissipative, and inertial, suggesting configuration d. In fact, it is by inertial coupling that portable phones, gaming pads or more advanced transducers operate [START_REF] Yao | Design and analysis of a recoil-type vibrotactile transducer[END_REF].

Option e can achieve an effect similar to that of arrangement c but introduces another tradeoff if the displacement of the distal motor is limited, since inertial coupling establishes a quadratic dependency between this displacement and the lowest frequency at which is can operate.

Dissipative forces can also be used to achieve desirable couplings. The counterpart of case a, case f, is described by [START_REF] Chew | Series damper actuator: a novel force/torque control actuator[END_REF] for robotic applications. It is not appealing for haptics since coupled to elastic load like a hand, this coupling will introduce an attenuation of 6 dB/octave in addition to the 12 dB/octave attenuation due to the motor inertia.

We could then think of yet another option, g, to 'take up' the response in the high-frequency, but similarly to b, this configuration would solve this problem by creating another. Option h is worth mentioning since it can be used to generate distortion-free haptic information up to 100 Hz but does not operate at dc [START_REF] Campion | Passive viscous haptic textures[END_REF]. Arrangement i would suffer from the same difficulties as e.

This enumeration naturally leads us to the last design, j, that could compete with the elastic coupling of configuration c. A comparison between the two options is carried out in the next section.

Discussion

Elastic couplings have the property of storing energy. This property is an advantage in applications where energy storage and restitution is a desirable property, such as in a walking machine, but in haptics, energy storage is inconvenient, even hazardous.

An elastically coupled drive system must be woundup to output a torque. If the output torque needs to be brought to zero abruptly, for instance when breaking contact with a virtual wall or when letting go of the handle during interaction, the limited speed of coupler unwinding forces a delay in the regulation of the output torque since this torque is two integrals away from the displacement of the motor. Moreover, the stored elastic energy could be uncontrollably converted to output kinetic energy, causing a hazard.

With viscous coupling the proximal stage must spin, also storing energy but in kinetic form. This energy cannot be transferred at the output transiently since the coupling naturally dissipates energy in itself, rather than in the user, without any need for active control. In contrast with the elastic coupling case, the energy storing quantity-velocity-is only one integral away from the output torque.

The basic servo mechanism that must govern the proximal stage is velocity control, which is easier to synthesize than a position servo since the system has one fewer order. The proximal stage can react more quickly than in the case of an elastic coupling, reducing the reliance on the distal stage to produce an accurate output torque.

There are other advantages of the viscous coupling over the elastic coupling, chiefly among them is accuracy. Dual-stage designs produce accurate output torques because, in essence, they are force-feedback systems. In the elastically coupled case, the feedback control system must regulate the spring deflection. If the elastic element is not of stringent quality, there could be hysteresis and high-frequency modes. More generally speaking, viscous coupling enjoys the general advantage of precluding the occurrence of structural dynamics present in any mechanical system having elastic elements.

With viscous coupling the onus is also on the coupler to provide an accurate relationship between dis-placement and torque. Due to the absence of contact, the main source of the injurious properties of mechanical systems, couplers based on eddy-current induction create resistive torques that are exactly proportional to the relative slip velocity. At very high velocities, the viscous coefficient decays, see the discussion in Section 3.2.1.

Prototype Drive

We now describe the system depicted in Fig. 2. The objective was to prototype a haptic interface that is such that, seen from the handle, the system's apparent mechanical impedance is small, but the torque that it can deliver is large. The result is a very highperformance system that closely meets the goals listed in Table 1 and thus is close to achieve perfect transparency.

Design Parameters

The quantities that participate in the performance are:

J 1 : the inertia connected to the proximal shaft that includes the proximal rotor and the clutch inductor;

J 2 : the inertia related to the distal axis that includes the rotor of the distal motor, the clutch armature, and the handle;

b: the coefficient of viscosity of the clutch.

Rapid changes in torque benefit from a short time constant for the proximal stage. This stage can be represented by a first-order system with a time constant equal to the ratio of its inertia to the viscosity of the clutch. Reducing inertia and increasing viscosity improves the capacity of the proximal stage to effect rapid torque changes without the help of the distal stage.

Achieving precise output torque requires the perturbing torque caused by friction to be minimized. In open loop, the inertias of both stages should be minimized. In closed loop, however, the inertia of the proximal axis can be hidden, since the energy needed to back-drive the system is supplied by an external source-the electronic amplifier, and not by the load.

The control is based on the estimation of the relative speed of the two shafts and can suffer from quantization effects. For a given torque transmitted by the clutch and a given quantization resolution, a higher viscosity decreases relative speed and thus the signal/noise ratio. These conflicting requirements require the use of high-resolution encoders or of a state observer.

In all cases, a design generally benefits from minimizing the inertia of the distal stage, J 2 , as well as minimizing the time constant of the proximal stage, J 1 /b.

Analysis

Design Criteria Approximations

The quantities J 2 and J 1 /b are affected by the mechanical and magnetic characteristics of the coupler. An eddy-current coupler has more parameters to optimize than a brake. There are two rotors: a permanent-magnet inductor and a disc armature. A multiple-loop magnetic circuit was designed according to the guidelines found in [START_REF] Gosline | Eddy current brakes for haptic interfaces: Design, identification, and control[END_REF] and from the results of the prototype described in [START_REF] Millet | The ultimate haptic device: First step[END_REF]. Figures 3a and3b illustrate the adopted geometry. As seen in Fig. 3b, a multiplicity of closed magnetic circuits create a periodic magnetic flux that crosses the annular armature mostly orthogonally. As a result, relative movement induces eddy-current loops as represented in Fig. 3c.

The analysis of the two optimization criteria, J 2 and J 1 /b, provided insight regarding the critical design parameters. The geometric parameters are defined in Fig. 3b and3c. The first criterion, J 2 , is the sum of the inertias of the armature, the rotor of the distal motor, and the handle. Provided that the width of the armature w a is much smaller than its mean radius r, to good approximation,

J 2 ≈ 2πρ a t a w a r 3 + J r2 , (1)
where ρ a is the mass density of the armature, t a is its thickness, and J r2 is the inertia of the distal rotor plus that of the handle. The mean radius of the clutch cubically scales the inertia of the armature, so it should be kept small. In our prototype, the handle accounts for about two thirds of J 2 and the armature for one third. The second criteria to minimize is the time constant of the proximal stage J 1 /b where 2n magnets, two iron annuli, the rotor of the proximal motor, and the magnet holder participate in the moment of inertia. Considering that the width of the annuli, w f , and the length of the magnets, l m , are much smaller than the mean radius, r, the inertia is approximated by

J 1 ≈ 2nl m w m t m r 2 ρ m + 4πw f t f r 3 ρ f + J r1 (2)
where ρ m and ρ f represent the mass densities of the magnets and of the annuli, t m and t f their thicknesses, respectively, and J r1 is the inertia of the proximal rotor plus that of the magnet support. The modeling of the magnetic field and induction in complex geometries is difficult because it involves several nonlinear phenomena [START_REF] Balakrishnan | Air-gap reluctance and inductance calculations for magnetic circuits using a schwarz-christoffel transformation[END_REF][START_REF] Heald | Magnetic braking: Improved theory[END_REF][START_REF] Srivastava | An alternative approach for calculation of braking force of an eddy-current brake[END_REF]. Nevertheless, several simplifying assumptions can be made to find a closed-form expression that estimates the viscosity of the clutch. When the armature rotates relatively to the magnets, the magnetic field, B, induces an electrical field, E = v × B, orthogonal to B and to the linear velocity v. Following the reasoning of [START_REF] Wouterse | Critical torque and speed of eddy current brake with widely separated soft iron poles[END_REF], an approach to estimate the current density induced in the armature is to suppose that the return path around the magnet is of negligible resistance. The current density induced under magnets is then J = σ a v × B, where σ a is the armature's conductivity. Assuming that B is uniform through the volume of the armature covered by magnets and is sufficiently greater than the magnetic field induced by the eddy currents, the power dissipated by the eddy currents can be computed by integrating the currents over the volume of the armature that is covered by magnets,

P diss = 1 σ a |J | 2 dV ≈ nσ a t a l m w m r 2 B 2 m ω 2 , ( 3
)
where n is the number of magnetic loops and B m the magnetic field produced by the magnets through the armature. The coefficient of viscosity, b, is then approximated by b ≈ nσ a t a l m w m r 2 B 2 m .

(4)

That the viscosity coefficient is independent of the velocity follows from the assumption than the magnetic field induced by eddy current can be neglected compared to B m , which is borne out at low speeds [START_REF] Anwar | A parametric model of an eddy current electric machine for automotive braking applications[END_REF]. For the present prototype, this assumption breaks down at high torques, that is, at speeds higher that 2 000 rpm where b would probably noticeably begin to decrease. This effect, however, is of no consequence to our main objective. We can now evaluate the proximal stage time constant

J 1 b ≈ 2nl m w m t m r 2 ρ m + 4πl m t f r 3 ρ f + J r1 nσ a l m w m t a r 2 B 2 m . (5)
From this expression, it follows that stronger magnets, higher conductivity, and a lower mass density in the armature all enhance performance, which is obvious. Finding the optimal number of magnets and their geometry, however, is not trivial.

Similar issues have been discussed for the design of magnetic undulators used in free-electron lasers and synchrotron radiation facilities. Literature from this topic provides some insight [START_REF] Quimby | High field strength wedged-pole hybrid undulator[END_REF][START_REF] Elleaume | Design considerations for a 1 Å SASE undulator[END_REF], but the tradeoff with inertia is not considered. These considerations justify a multi-physics numerical design optimization approach.

Numerical Simulation and Optimization

A 3D simulation was set up with a multiphysics simulation package (comsol ab, Stockholm, Sweden), a finite element analysis simulator suitable to simulate eddy currents induced in a moving conductor.

The two symmetries (mid-plane and n angular antiperiodicity) in the problem were used to reduce computation time. By establishing anti-periodic boundary conditions at the mid-planes of two consecutive magnets, the simulation of one of the 2n anti-periodic sectors was sufficient. The mesh resolution was set to 0.3 mm for the regions of the armature shadowed by the magnets and coarser elsewhere. An important criterion for selecting the mesh resolution was to mesh the boundary of the annulus with element sizes sufficiently small to account for the skin effect. An ungauged potential formulation was used to solve the electromagnetic problem. Typically, convergence was achieved in 15 steps for a relative precision of 10 -6 .

The skin depth was about 12 mm for a relative speed of 21 rad/s, hence invading the whole armature at higher speeds. This effect is visible in Figure 4 that shows the current density induced in the armature. Density is slightly smaller in the mid-plane of the armature. Figure 4c shows the results of the simulation for a relative speed of 21 rad/s, generating a braking torque of 50 mN•m. The influence of the thickness and width of the armature, t a and w a , the number of pairs of magnets n, and the gap, g, between the armature and magnets were studied. Other parameters were fixed. The mean radius was r = 20 mm, the magnet dimensions 2×2×10 mm, their remanent field B m = 1.37 T. The geometry of the iron annuli was constrained by the magnets and the need to prevent magnetic saturation.

Figure 5 shows that increasing the number of magnets increases the total torque until a certain limit where magnets become too close, causing magnetic leakage. The reduction of the inertia of the proximal stage J 1 conflicts with the optimal number of mag-nets. The best tradeoff for the present design corresponds to a number of magnets pairs between 32 and 36 and an armature thickness of 1.0-1.5 mm. Figure 6 shows that selecting the width of the armature to be slightly longer than the length of magnets significantly increase the viscosity coefficient, while not affecting the inertia of the distal axis significantly. This result agrees with experimental observations found in [START_REF] Gosline | Eddy current brakes for haptic interfaces: Design, identification, and control[END_REF]. Lastly, Figure 6 shows that the viscosity coefficient is strongly influenced by the gap which should be kept as small as possible, within the limits of manufacturing feasibility.

Resulting Prototype

The system seen in Figure 2 was manufactured using a combination of conventional machining and rapid prototyping (abs plastic), which was effective for those parts which had to be lightweight. Some structural engineering was carried out to ensure that all the parts had the required rigidity. The motors were a 35 mm core-less dc motor (model re-35, 273754, graphite brushes, Maxon Motors AG, Sachseln, Switzerland) for the proximal stage and a small 16 mm motor (model re-16, 118698, precious metal brushes) with a peak torque approximatively 20 times smaller than the larger motor. The motors were powered by two analog current/voltage amplifiers (model lcam, Quanser, Markham, on, Canada). The two shafts were fitted with high-resolution digital encoders (Mercury M1800, MicroE Systems, Bedford, ma, usa) giving 1,638,400 counts per turn. The clutch had 60 neodymium magnets, and the armature parameters were t a = 1 mm, r = 20 mm, w a = 12 mm.

It is of paramount importance to have accurate estimates of the system parameters. They were measured as described in [START_REF] Millet | The ultimate haptic device: First step[END_REF] and the results collected in Table 2. The viscosity coefficient was a bit lower than expected from the simulation because of the lack of dimensional accuracy of the rapid prototyping manufacturing process forcing us to increase the gap.

-3 N•m/A viscous coefficient b 8.7 10 -4 N•m•s/rad proximal inertia J 1 2.6 10 -5 kg•m 2 distal inertia J 2 6.4 10 -6 kg•m 2 distal dry friction C 2fric 1.3 10 -4 N•m max torque C 1max 2.0 10 -1 N•m
With proper control, the system matched closely the entire range of human sensorimotor performance, or "perfect transparency" as quantified in Section 2.

Control System Design

Referring to Fig. 7, the system has two torque commands and one torque output, which can potentially provides for a number of control options. Here, the primary control objective is output torque regulation and tracking. One effect of this particular control objective is to disconnect the effects of the proximal motor dynamics from the experience of the user. During active exploration, when the handle is back-driven, the power required to move the proximal motor is entirely supplied by the power amplifier, not by the user.

To achieve this objective, the reference torque, τ d , is compared to the torque produced by the coupler, which is proportional to relative velocity between the proximal and distal motors. The proximal motor is slaved by the compensator C(z) to the torque error, operating in the sampled-data domain with period h. The distal motor is commanded to reduce the error without compensation since it is needed to transiently compensate for the slower response of the proximal stage. More complex schemes to achieve a wider range

z -1 z -1 h h ZOH ZOH C(z) τ d τ user 1 s 1 s b -b θ1 b h τ out τ coupler θ1 θ2 1 b + J 2 s 1 b + J 1 s θ2 θ1
Figure 7: Control Scheme.

of objectives could involve compensation of the distal motor dynamics as well.

The controller has a path to regulate the differential velocity of the two shafts ( θ1 -θ2 ). If regulation and tracking are sufficiently good, the namics of the proximal stage is from the user's experience. The transparency is limited by the high-frequency dynamics of the distal stage only. In other words, for a null torque set-point, what is felt is only the distal stage. This scheme relies on an accurately specified viscous coefficient, b, which eddycurrent couplers allow.

Polynomial Pole-Placement

We adopted a discrete-time polynomial pole placement approach because, in contrast with other approaches, this method allows for an exact design [START_REF] Landau | Differential elastic actuator for robotic interaction tasks[END_REF][START_REF] Åström | Computercontrolled systems: theory and design[END_REF][START_REF] Ostertag | RST-controller design for sinewave references by means of an auxiliary diophantine equation[END_REF]. The discrete-time controllers are specified by three polynomials R(z), S(z), and T (z), as shown in Fig. 8, where y(z) is the measured plant output, r(z) the reference set-point, and u(z) the control signal. Control synthesis is purely algebraic and brings about a solution systematically, once the discrete model of the plant is known and a closed-loop model transfer function, H m (z), has been specified. The three control polynomials can be selected to achieve certain performance objectives, including regulation and tracking. The discrete-time plant, here the sam-pled continuous-time plant, is

H(z) 1/S(z) T (z) R(z) u(z) r(z) y(z) H m
H(z) = B(z) A(z) , (6)
and the control law is

U (z) = T (z) S(z) r(z) - R(z) S(z) y(z). (7)
The feedback control, -R(z)/S(z), gives a closedloop characteristic equation,

A(z)S(z) + B(z)R(z) = P (z), (8)
that takes the form of a Diophantine equation. The design involves first selecting P (z) to determine the closed-loop regulation behavior through the placement of the closed-loop poles, taking the constraints on the degrees of the polynomials into account. The tracking polynomial, T (z), is then selected to independently determine the tracking behavior.

Design

If we neglect the velocity of the distal stage, that is, consider it small compared to that of the proximal stage, then the compensator acts on the dynamics of the proximal stage only. The discrete-time transfer function of the plant results from cascading a zero-order hold with an integrator and discrete-time derivation. With a sample period of 0.1 ms, the plant transfer function becomes,

H(z) = 0.001671z + 0.001669 z 2 -0.9967z . (9)
The desired closed-loop polynomial was selected to correspond to a second order continuous plant with natural frequency and damping, ω 0 = 680 rad/s, ζ = 0.9, respectively. This led to the characteristic polynomial

P (z) = z 2 -1.8800z + 0.8848. ( 10
)
The solution of the Diophantine equation ( 8) gave

S(z) = z -0.98855, R(z) = 1.3240. ( 11
)
Perfect model tracking is achieved when

T (z) = P (z) B(1) , (12)
and the model was

H m (z) = 0.001625z + 0.001569 z 2 -1.898z + 0.9009 . ( 13
)

Control Performance

The control system was implemented on a personal computer running rtai real-time Linux kernel and the sampling rate was set to 10 kHz. The system was tested by holding the handle with a normal grip and demanding a step of torque. Saturation levels were set to 5 mN•m for the distal motor and to 200 mN•m for the proximal motor, because the distal stage motor normally operates only during the transients, its saturation can be set much higher, viz. by a factor five, without adverse effects. Figure 9a shows the response of the system to a step of amplitude 5 mN•m, where the output torque, τ out (black line), is the sum of the viscous coupler torque, τ coupler , and of the distal stage motor torque, τ 2 . As was expected, see Fig. 9a, the proximal motor received a large command, τ 1 , during a short transient to minimize the reaction time and the torque missing in τ coupler was taken up by the torque, τ 2 , produced by the distal stage motor. When saturation was not reached for the distal stage motor, as in the case of Fig. 9a, the mechanical time constant of the system was equal to that of the distal motor, that is about, 1.0 ms. In the cases where the demanded torque caused saturation of the distal stage, the transient response was imperfect, as in Fig. 9b where the step magnitude is twice that of the distal motor saturation level. The steady-state part, however, remains unaffected.

Validation

The objective was to verify that the high level degree of transparency of the dual-stage drive had a measurable effect on human performance when using a haptic interface. The purpose of the experiment was to test the device effectiveness, not human performance. To this end, we selected a task where participants had to perform as signal detectors using a stimulus that had two components. The first component had a fixed magnitude and the other was varied from imperceptible values to easily detectable values. To create an experimental condition that resembled actual use, we asked the same volunteers to detect the presence of weak, high-frequency force perturbations superposed onto a slowly varying background force at random locations in the workspace. At any given time, however, the participants were not aware of the condition under which they were working.

The stimuli were administered through two identical handles, one connected to the proximal motor, and other to the distal motor, see Fig. 10. We tested whether the subjects could detect smaller details in the same signal rendered with the dual-stage interface (Condition A) rather than with the large motor alone (Condition B), similar to a conventional interface. All other parameters remained equal.

Figure 10: Experimental setup. In Condition A, the handle was connected to the output of the dual-stage drive. In Condition B, an identical handle was connected directly to the shaft of the proximal motor and the coupler was disconnected. The device was mounted on a heavy steel slab to prevent possible vibrations arising from the rotation of the proximal motor to find a path to the hand.

Methods

Subjects

Six members of the laboratory, five male and one female, volunteered their time. They were aged from 23 to 56 with a median age of 28.

Stimulus

The slowly varying signal was produced by simulating dry friction and the fast signal was obtained by modulating the amplitude of the force according to

   τ C 1 -h 0 sin 2π θ λ , when θ 0 < θ < θ 0 + Θ, τ C , otherwise, (14)
see Figure 11. The friction torque was calculated according to the algorithm described in [START_REF] Hayward | A new computational model of friction applied to haptic rendering[END_REF] and its level, τ C , was fixed at 5 mN•m (0.14 N at the handle periphery). From trial to trial, the onset of the texture, θ 0 , changed randomly in the range 0-1.57 radian (0-π/2), and its extent, Θ, was fixed at 0.035 radian (1.22 mm) and its period, λ, at 0.005 radian (0.17 mm). When turning the handle, it felt as if its motion was impeded by a brake and it was like trying to detect a slight imperfection. For Conditions A and B, preliminary tests enabled us to determine that the thresholds of detection for h 0 were below 0.1 and 1.0, respectively (0.014 N and 0.14 N at the handle periphery). The period, λ, and amplitudes, h 0 , simulated very fine and very small textures, which was a test for the device's capabilities. A better device should allow users to perform better.

Rendering a given texture with a given device, however, imposes limits that are not to be exceeded [START_REF] Campion | Fundamental limits in the rendering of virtual haptic textures[END_REF]. The lowest spatial period, λ, that can be reliably synthesized is constrained by (α/λ)vT < 1 where v is the velocity and α > 2, and T is the sampling period. For a velocity of at most 2.9 rad/s (0.1 m/s at the handle periphery), with T = 10 -4 s and α = 10, the smallest spatial period should be 2.9 mrad (0.1 mm). The other constraints in were easily met by our prototype. The width, Θ, of the texture should also be greater than 3λ.

Procedure

We adopted a two-alternative forced choice testing procedure where h 0 was randomly selected in the sets {0.0, 0.01, 0.03, 0.053, 0.08, 0.1} mN•m and {0.0, 0.2, 0.4, 0.6, 0.8, 1.0} mN•m, for Condition A and B, respectively. Each stimulus was presented 10 times, therefore the participants completed 60 trials for each condition, which could be done in a few minutes. In both conditions, the setup was concealed inside an enclosure having the same visual aspect from both sides. Only the two handles were visible. To switch from Condition A to Condition B, all there was to do was to turn the set-up around and to disconnect the coupler from the proximal motor.

The participants were asked if they felt anything at all, in which case they had to respond 'no', or a small texture, in which case they would answer 'yes', and indicate where they had felt it.

Results

We fit the experimental data with cumulative Gaussian distributions of the form f (x) = [1 + erf((h 0µ)/ √ 2σ)] where µ and σ 2 were the detection mean and variance, respectively, using a non-linear leastsquare procedure. The threshold points were extracted from the amplitude of the stimulus texture that corresponded to a 50% probability texture detection, for each subject, in the configurations. The average detection threshold across subjects was 0.036 when using the dual-stage drive and 0.27 when using the standard configuration.

Discussion

All by-and-large, produced typical psychometric curves shown in Fig. 12. Participant SD, who was very cautious when using the improved device, made decisions without hesitation when using the standard drive, but the opposite was true for participant SE. Such differences are not surprising considering that haptic interfaces, no matter how well they are constructed, are instruments and that users can choose to focus on various aspects of their operation, even in a simple detection task. They can also refer to their past experience to various degrees. Whether naive or experienced, their relationship with the instrument and their skill in using it can vary and the effects of changing an instrument design can be different from person to person. Despite these inter-subject variations, it can then be concluded that the performance detection of the subjects was improved by almost an order of magnitude. It was the case for participants SA, SB, SC, and SF) and to a slightly lesser extend for SD and SE. Lastly, the thresholds corresponded to force variations of 5.2 mN and 38.3 mN for Conditions A and B respectively, the former figure being in line with the human performance threshold (0.5 g experienced dynamically).

Conclusion

Using a dual-stage approach, we designed an electromagnetic drive for use in haptic interfaces and other applications where the maximum torque delivered is decoupled from its effective inertia. In terms of haptic interfacing, this means that transparency is improved since the 'tool' interposed between the source of torque and its point of delivery causes parasitic forces originating from inertia and friction to be commensurate with perceptual thresholds. In other words, the torque specified is for all practical purposes independent from both the movements of the hand and from those of the device.

In sum, the interface achieves two orders of magnitude of improvement in transparency over existing designs measured in terms of the magnitude of parasitic forces owing to friction and inertia. The inertial load seen by a finger is only of the order of one or two grams. Yet, the design does not limit the torque magnitudes that can be produced, and can operate in its present configuration over a bandwidth that exceeds one kilohertz. For the task that we have selected, the improved transparency allowed users to detect details that were ten times smaller in magnitude than when using a conventional design.

Figure 1 :

 1 Figure1: Several ways one or two motors can been coupled to a load. See text for discussion. Combination j is adopted for the device described in the present article.

Figure 2 :

 2 Figure 2: Prototype realized according to strategy j.

Figure 3 :

 3 Figure 3: Eddy-current viscous coupler. a, Cutaway view of the coupler. b, Magnetic flux density streamlines generated by the magnets. c, Current density streamlines in the mid-plane of the armature during relative movement.

Figure 4 :

 4 Figure 4: a, Current density induced in a 32 magnet pair armature rotating at 21 rad/s. b, Zoom on one sector. c, Electromagnetic simulation of a half magnetic loop.

Figure 5 :

 5 Figure 5: Viscosity coefficient and time constant of the proximal stage axis according to the armature thickness and the number of magnets. Other parameters were fixed, w a = 10 mm and g = 1 mm.

Figure 6 :

 6 Figure 6: Influence of the gap and the width of the armature on the viscosity coefficient.

Figure 8 :

 8 Figure 8: Polynomial control structure.

Figure 9 :

 9 Figure 9: Step response time domain performance. a, Case without saturation. b, with saturation: the distal motor torque, τ 2 , is limited to 5 mN•m.

Figure 11 :

 11 Figure 11: Stimulus signal. A small burst of oscillation of size Θ, amplitude h 0 , and spatial period λ, was superposed onto a larger friction force, τ C .

Figure 12 :

 12 Figure 12: Results and fits for Conditions A and B.

Table 2 :

 2 Measured system parameters.

	Quantity	Value
	distal torque constant k 2	7.9 10

for brevity, no distinction is made between forces and torques.

Acknowledgments

The authors wish to thank Andrew H. Gosline for advice on magnetic design of viscous couplers and Bruno Sauvet for help with prototype manufacturing. This work was supported by the ANR (Agence Nationale de la Recherche, France), project PACMAN "perception haptique des échelles micro et nanoscopiques". Additional funding was from the European Research Council, Advanced Grant PATCH, agreement No. 247300 (to V.H.).

Authors' submitted copy, published in the International Journal of Robotics Research, 31(3):319-329, 2012.