D.M.D Xiaoqian Ye

Ph.D Audrey Guilmatre

Ph.D Boris Reva

Ph.D Inga Peter

Ph.D Yann Heuzé

Ph.D Joan T Richtsmeier

Deborah J Fox

Rhinda J Goedken

M.D Ethylin Wang Jabs
email: ethylin.jabs@mssm.edu

Ph.D Paul A Romitti
email: paul-romitti@uiowa.edu

Mutation Screening of Candidate Genes in Patients with Nonsyndromic Sagittal Craniosynostosis

Background-Craniosynostosis is a condition that includes the premature fusion of one or multiple cranial sutures. Among various craniosynostosis forms, midline sagittal nonsyndromic craniosynostosis (sNSC) is the most prevalent. Although different gene mutations have been identified in some craniosynostosis syndromes, the etiology of sNSC remains largely unknown.

INTRODUCTION

Craniosynostosis is defined by the premature fusion of one or more cranial sutures and the malformation of cranial bones that result in abnormal skull shape. It is one of the most common birth defects, occurring in 1 per 2,000 to 2,500 live births. Craniosynostosis patients are often classified according to the suture(s) involved and whether the cranial vault features occur as an isolated defect (nonsyndromic form), or are associated with other facial, specific limb, or dermatological features as part of a syndrome (syndromic form). Nonsyndromic craniosynostosis (NSC) typically involves premature closure of a single suture, the most frequent type being sagittal, and it is controversial whether it is followed by the coronal (unicoronal being more frequent than bicoronal synostosis) or the metopic suture. [START_REF] Heuzé | Closing the gap: Genetic and genomic continuum from syndromic to nonsyndromic craniosynostoses[END_REF] The premature fusion of the sagittal suture results in an increase of the anteroposterior direction of the skull, particularly the calvarium, with growth restriction in the transverse direction, causing a narrow, elongated, boat-shaped skull deformity (dolichocephaly or scaphocephaly) and varying degrees of frontal bossing and occipital bossing. Abnormal growth of the skull can subsequently lead to increased intracranial pressure that can be associated with changes in brain morphology, which in turn can contribute to visual impairment and neurocognitive deficiency. [START_REF] Zeiger | Genetic and environmental risk factors for sagittal craniosynostosis[END_REF][START_REF] Boyadjiev | International Craniosynostosis Consortium. Genetic analysis of non-syndromic craniosynostosis[END_REF] The treatment of sagittal craniosynostosis is currently evolving. Total vault remodeling (CVR) has been effective for correcting cranial shape and expanding intracranial volume. Minimally invasive techniques, such as endoscopic suture release, spring assisted surgery, and distraction osteogenesis have been used in an effort to potentially reduce surgical morbidity. [START_REF] Gerety | Operative management of non-syndromic sagittal craniosynostosis: A head-to-head meta-analysis of outcomes comparing three techniques[END_REF] Sagittal nonsyndromic craniosynostosis (sNSC) accounts for approximately one-half of all affected patients and shows a strong male predominance (male:female ratio of 3.5:1). [START_REF] Lajeunie | Genetic study of scaphocephaly[END_REF] Segregation analysis indicates that sNSC follows an autosomal dominant inheritance with a reduced penetrance in only 6% of cases. [START_REF] Lajeunie | Genetic study of scaphocephaly[END_REF] The majority of sNSC is thought to be sporadic. The current genetic understanding suggests that the underlying etiology of sNSC is complex, most likely involving both genetic and environmental exposures (broadly defined). [START_REF] Kallen | Maternal smoking and craniosynostosis[END_REF] In the past two decades, intensive investigations among the more than 180 craniosynostosis syndromes have provided valuable insights into identifying genes involved in calvarial suture development. [START_REF] Eacute | Genetics of craniosynostosis: Genes, syndromes, mutations and genotype-phenotype correlations[END_REF] Among the identified syndromic-causing genes, some are responsible for rare syndromes with midline sagittal or multiple suture craniosynostosis (Table 1). Heterozygous gain-of-function mutation hotspots in fibroblast growth factor receptors 1, 2 and 3 (FGFR1, FGFR2, FGFR3) account for most of the common craniosynostosis syndromes, including Apert, Crouzon, Pfeiffer, Jackson-Weiss, Beare-Stevenson and Muenke syndrome. [START_REF] Melville | Genetic basis of potential therapeutic strategies for craniosynostosis[END_REF] In addition to those FGFR-related craniosynostosis syndromes, mutations in several transcription factors such as TWIST1 and MSX2 have been implicated in Saethre-Chotzen syndrome and Boston-type craniosynostosis respectively. The clinical phenotypes of these syndromes vary greatly, but all have coronal suture involvement. In contrast to the increasing number of human genes associated with syndromic craniosynostosis, [START_REF] Melville | Genetic basis of potential therapeutic strategies for craniosynostosis[END_REF] only variants in a small number of genes, FGFR1-3, TWIST1, ALX4, ERF, IGF1R, and LRIT3 (Table 2) have been reported in sNSC patients. The mutation hotspots that are important for syndromic craniosynostosis, FGFR1 exon IIIa, FGFR2 exons IIIa and IIIc, FGFR3 exon IIIa, and TWIST1 exon 1 and their somatic changes, have been less likely to be associated with sNSC, 2,3,9,10 perhaps because the sagittal suture is infrequently closed prematurely in syndromic patients. As such, patients with nonsyndromic unilateral or bilateral coronal synostosis merit the genetic testing of mutation hotspots in FGFR1-3 and TWIST1, but mutation testing is not proposed for patients with sNSC unless a syndromic association is suspected due to familial recurrence. [START_REF] Boyadjiev | International Craniosynostosis Consortium. Genetic analysis of non-syndromic craniosynostosis[END_REF][START_REF] Roscioli | Genotype and clinical care correlations in craniosynostosis: findings from a cohort of 630 Australian and New Zealand patients[END_REF] With regard to sNSC, a genome-wide association study (GWAS) showed significant association in a 120kb region downstream of BMP2 and within a 167kb region of BBS9 in 130 non-Hispanic white caseparent trios (N=130) and replicated in 172 case-controls of the same ethnicity, but without identifying plausible causative mutations. [START_REF] Justice | A genome-wide association study identifies susceptibility loci for nonsyndromic sagittal craniosynostosis near BMP2 and within BBS9[END_REF] To further investigate the pathogenesis of sNSC, we sequenced known FGFR1-3 and TWIST1 mutations that were found to be related to sNSC (Table 2) and two new candidate genes RAB23 and BMP2 using case-parent trio specimens from families identified from two population-based birth defect surveillance programs in the U.S.

MATERIALS AND METHODS

Study Population and Data Collection

Cases were live born children with a diagnosis of single suture, sagittal craniosynostosis born from January 1, 1998 through December 31, 2008 to resident mothers in Iowa or New York State. Cases were enumerated from the Iowa Registry for Congenital and Inherited Disorders and the New York State Congenital Malformations Registry, two populationbased birth defect surveillance programs. Controls were live children born without a major birth defect diagnosis to resident mothers in Iowa or New York State during the same time frame as cases and frequency matched by birth year to cases. To be eligible for recruitment, a case or control mother needed to be able to complete an English language questionnaire and to have custody of the case or control child.

Data collection was conducted in three phases: family history questionnaire, saliva specimen collection, and medical record abstraction. A case or control mother was sent a pre-contact letter informing them about the study, and two weeks later, an introductory packet comprised of an introductory letter, the family history questionnaire, a fact sheet/rights as a research subject, a frequently asked questions sheet, and $10 reimbursement was mailed. A systematic follow-up protocol was used to encourage return of a completed family history questionnaire. Following receipt of a completed questionnaire, a case or control mother was sent a saliva specimen collection kit comprised of: a specimen collection letter and instructions; written consent forms for mother and father; specimen collection supplies for mother, father and child; $20 reimbursement; and a medical release form (case mothers only). Again, a systematic follow-up protocol was used to encourage return of saliva specimens and for case mothers, a medical release form. Once a family returned the saliva specimens, a thank you letter and $20 were sent to the mother. For case mothers who returned a signed medical record release form, medical record abstraction for postnatal events, including request for computed tomographic cranial images, was completed. The study protocol was approved by the Icahn School of Medicine at Mount Sinai, The University of Iowa and New York State Department of Health institutional review boards, and conformed to the guidelines set forth by the Declaration of Helsinki.

Mutation Screening

Candidate genes or mutations were selected based on their known involvement in human sagittal craniosynostosis, in animal models, and/or expression studies related to sagittal suture development. In addition to selected variants in FGFR1 (c.1508C>T), FGFR2 (c.833 G>A; c.943 G>A), FGFR3 (c.749C>G; c.1000 G>A) (Table 2), and the TWIST1 gene, we also included in our screen two new candidate genes: RAB23, which is responsible for Carpenter syndrome, which presents with midline sagittal, as well as metopic, synostosis, and BMP2, which was identified from the previously mentioned GWAS study. [START_REF] Justice | A genome-wide association study identifies susceptibility loci for nonsyndromic sagittal craniosynostosis near BMP2 and within BBS9[END_REF] DNA sequence analysis was performed on 93 affected sNSC cases. Primers for the coding regions of all candidate genes or variants (FGFR1 NM_023105; FGFR2 NM_000141; FGFR3 NM_000142.2; TWIST1 NM_000474; RAB23 NM_183227; BMP2 NM_001200) were designed using Primer 3 (http://bioinfo.ut.ee/primer3-0.4.0/) (see Supplemental Digital Content 1, Table shows sequencing primers for candidate genes, INSERT LINK). The sequencing was conducted by Genewiz Inc (http://www.genewiz.com/) and sequences analyzed using Sequencher software (v4.8; Gene Codes Corporation, Ann Arbor, MI). We initially screened only the case child. When a variant was identified, we sequenced maternal DNA specimens, and if available, paternal DNA specimens to determine if the variant was de novo or segregated within the family. To confirm the variation, we carried out molecular analysis of 58 healthy children, or 116 control alleles, and compared our results to the 1000 Genomes Project (1KG) database (http://browser.1000genomes.org), the NHLBI Grand Opportunity Exome Sequencing Project (ESP) database (http://evs.gs.washington.edu/ EVS/), the Exome Aggregation Consortium (ExAC) database (http:// exac.broadinstitute.org/) and Locus Specific Mutation Database (LSDB) (http:// grenada.lumc.nl/LSDB_list/lsdbs).

RESULTS

Among the 93 affected sNSC cases screened, we identified two rare novel coding variants in our case child sample. These variants included a heterozygous c.439C>G transversion in TWIST1 exon 1 (Fig. 1), which predicted a glutamine-to-glutamic acid transition that occurred within the highly conserved loop domain. This variant was also present in the case child's mother, who had an undefined jaw anomaly with corrective surgery when she was a teenager. The other heterozygous variant, c.730_731insG (NM_023105), was observed in exon 6 in FGFR1 in a different case child (Fig. 1), but not in the case child's mother; a specimen was not available for the case child's father to be able to determine if this variant was a de novo variant. This insertion was predicted to abolish the entire immunoglobin III domain, including the ligand binding region. These two variants were not seen in the 1KG, ESP or LSDB databases, however the TWIST1 variant c.439C>G was found in an extremely low frequency in ExAC database (1 in 120318 alleles, phenotypic information unavailable). In addition, we identified a heterozygous variant, c.546A>C, in RAB23 in a third case child; this variant predicted the amino acid substitution of glutamic acid to aspartic acid. There is not sufficient extant data in the literature to predict a functional effect of this variant from evolutionary conservation using Mutation Assessor tool. [START_REF] Reva | Predicting the functional impact of protein mutations: Application to cancer genomics[END_REF] The case child's unaffected mother carried the same variant with no obvious phenotype, and this rare variant was also found in one individual from the NHLBI ESP database. It is unknown if the carrier from ESP had any relevant phenotype. Also we did not find any of the previously reported variants in FGFR1-3, TWIST1, RAB23 and BMP2 in our study samples (Table 2).

DISCUSSION

Using a population-based sample, we identified a very rare variant glutamine-to-glutamic acid transition in exon 1 in TWIST1 (c.439C>G; Fig. 1). The assessment of the functionality of this variant made by the Mutation Assessor tool is presented in Figure 2 and Figure 3. [START_REF] Reva | Predicting the functional impact of protein mutations: Application to cancer genomics[END_REF] From multiple sequence alignment, we concluded that the variant residue occupies one of the top specificity positions (Fig. 2) and contributes to the complex formation between transcriptional regulator TWIST1 and the transcriptional factor E2 alpha (Fig. 3). Interestingly, both types of the variant residues, the original glutamine and variant glutamic acid are present in specific subfamilies. This suggests that the glutamic acid variant may be responsible for the specific regulation of the DNA transcription, biologically different from the one performed by the common glutamine variant of the TWIST1. We believe this variant to be pathogenic based on the following evidence. First, this variant was located within the highly conserved loop domain, which is important for the TWIST1 protein's tertiary structure and dimerization. A number of mutations, including c.435G>C and c.443C>A, have been detected in different Saethre-Chotzen patients within the coding sequence for the loop domain at the surrounding nucleotide positions and were confirmed to be deleterious. [START_REF] Rose | The TWIST gene, although not disrupted in Saethre-Chotzen patients with apparently balanced translocations of 7p21, is mutated in familial and sporadic cases[END_REF][START_REF] Kress | Saethre-Chotzen syndrome caused by TWIST 1 gene mutations: Functional differentiation from Muenke coronal synostosis syndrome[END_REF] Second, the variant was extremely rare in the general population. It was absent in a panel of 116 normal control alleles and two of the largest human genome databases; and it was found only once in more than 6000 individuals in ExAC database and the individual information is unknown. The low frequency of the variant could not exclude its pathogenicity. In previous studies, some mutations which were found in other large cohorts of sNSC patients or confirmed with functional analysis also showed a very low population frequency in the public databases (Table 2). Lastly, the variant was found in both the case child and mother. Although the mother had no features of obvious craniosynostosis, she was reported to have anomalous jaw formation, another developmental bone abnormality. TWIST1, a helix-loop-helix transcription factor, has been implicated both in Saethre-Chotzen syndrome and in nonsyndromic coronal craniosynostosis. [START_REF] Roscioli | Genotype and clinical care correlations in craniosynostosis: findings from a cohort of 630 Australian and New Zealand patients[END_REF][START_REF] Wilkie | Prevalence and complications of single-gene and chromosomal disorders in craniosynostosis[END_REF] However, detected point mutations in TWIST1 have been implicated in very few cases of sNSC (Table 2) and are likely to contribute to very few patients. Variable phenotypic expressivity in patients with TWIST1 mutations is common. A S188L mutation in the TWIST Box domain was previously identified to be pathogenic both in a patient with sNSC and the patient's father, who had only small, square-shaped ears and is considered to represent a very mild manifestation of the TWIST1 mutation. [START_REF] Seto | Isolated sagittal and coronal craniosynostosis associated with TWIST box mutations[END_REF] Similarly in our study, the mother might have variable expressivity of the phenotype with only jaw defects. In the Twist1 +/-mice, sagittal suture closure was found to be related to the canonical Wnt signaling downregulation and the causative effects for sNSC are most likely of an epigenetic nature. [START_REF] Behr | Differential activation of canonical Wnt signaling determines cranial sutures fate: A novel mechanism for sagittal suture craniosynostosis[END_REF] Few FGFR1 mutations have been associated with craniosynostosis. A specific mutation P252R in Pfeiffer syndrome and a heterozygous FGFR1 I300L were found in a trigonocephaly patient with isolated metopic synostosis. [START_REF] Kress | An unusual FGFR1 mutation (fibroblast growth factor receptor 1 mutation) in a girl with non-syndromic trigonocephaly[END_REF] With our population-based sample, we also screened the previously reported FGFR1 c.1508C>T variant, [START_REF] Seto | Isolated sagittal and coronal craniosynostosis associated with TWIST box mutations[END_REF] but failed to detect the same variant in our study sample. Instead, we identified another novel variant, c. 730_731insG, near the c. 1508C>T variant, which led to a premature stop codon (p.A244fs*26) and was predicted to abolish the entire IgIII domain, including the ligand binding region. The function of c.730_731insG needs to be further studied. Gain-of-function mutations in FGFR2 have been identified in a majority of autosomal dominant craniosynostosis syndromes that usually include the coronal suture, but multiple affected sutures (coronal and sagittal) can also be involved in some complex cases. In addition to the mutation hotspots found in syndromic craniosynostosis, a few FGFR2 mutations have been reported in different patients with nonsyndromic sagittal craniosynostosis (Table 2). In contrast to FGFR2, only four mutations in FGFR3 have been identified in patients with craniosynostosis, P250R mutation in Muenke syndrome, A391E in Crouzon syndrome with acanthosis nigricans, P250L in one family with nonsyndromic craniosynostosis (sutures involved were unknown), and A334T in a family with mild sagittal crainosynostosis. [START_REF] Boyadjiev | International Craniosynostosis Consortium. Genetic analysis of non-syndromic craniosynostosis[END_REF][START_REF] Schindler | Heterozygous P250L mutation of fibroblast growth factor receptor 3 in a case of isolated craniosynostosis[END_REF][START_REF] Barroso | Mild isolated craniosynostosis due to a novel FGFR3 mutation, p.Ala334Thr[END_REF] A complete screening of FGFR3 was undertaken in 97 craniosynostosis patients with no known molecular defect and no causative mutation was identified. [START_REF] Wilkie | Prevalence and complications of single-gene and chromosomal disorders in craniosynostosis[END_REF] Our screening for FGFR2 C278Y and A315T variants and the FGFR3 mutation hotspots P250R and A334T failed to identify any novel sequence variants in either gene.

Carpenter syndrome is a rare craniosynostosis syndrome with midline sagittal suture involvement. Presumed loss-of-function mutations were found in RAB23, which encodes a member of the RAB guanosine triphosphatase family of vesicle transport proteins and acts as a negative regulator of hedgehog signaling. We screened the entire RAB23 and BMP2 coding and flanking regions in our study sample, but failed to identify any novel sequence variants in either gene. Our study is also the first population-based cohort that demonstrates that mutations in several genes associated with currently known or identified craniosynostosis usually associated with coronal craniosynostosis, were not common in our sNSC cases. Our sequencing results were consistent with the data from previous reports in non-population, Caucasian cohorts. [START_REF] Zeiger | Genetic and environmental risk factors for sagittal craniosynostosis[END_REF][START_REF] Boyadjiev | International Craniosynostosis Consortium. Genetic analysis of non-syndromic craniosynostosis[END_REF][START_REF] Justice | A genome-wide association study identifies susceptibility loci for nonsyndromic sagittal craniosynostosis near BMP2 and within BBS9[END_REF] Thus, sNSC remains a heterogeneous and challenging disorder. Furthering knowledge of gene variants that produce craniosynostosis will require identifying functional variants, investigating pathway and other interactions by using whole exome/genome sequencing, and conducting genome-wide association studies to identify genetic interactions that may reveal critical genetic relationships that modify craniosynostosis phenotypes.

CONCLUSIONS

The etiology of sNSC is complex and is yet unknown. Our sequencing results of a population-based cohort of sNSC patients indicate that the FGFR1-3 mutational hotspots, TWIST1, RAB23 and BMP2 gene coding regions which are involved in syndromic sagittal craniosynostosis and/or sagittal nonsyndromic craniosynostosis contribute to only a few cases of sNSC (variant detection frequency about 1%). Additional studies are needed to elucidate the pathological mechanisms underlying sNSC since these candidate genes may not be the major causes of sNSC and genetic testing of these genes was negative in most cases. Assessment of the functional impact of the variant p.Q147E (c.439 C>G) in TWIST1. A reduced presentation of the multiple sequence alignment of the protein region with mutated residue produced by Mutation Assessor. [START_REF] Reva | Predicting the functional impact of protein mutations: Application to cancer genomics[END_REF] The whole alignment is composed of 393 sequence homologs of TWIST1 (both human and non-human) and divided into 29 subfamilies; only 50 sequences and 4 subfamilies are shown. The positions in the alignment are assessed by evolutionary conservation within entire protein family and by specificity, i.e. evolutionary conservation within subfamilies; the lengths of the bars above the alignment represent the relative strength of conservation and specificity in positions of the multiple alignment. The residues, which are conserved within subfamilies, i.e. specific to the subfamilies, are typically responsible for the actual binding specificity of proteins. The mutated residue is located in one of the top specificity positions, both the original glutamine and the mutated glutamic acid residues are represented by subfamilies; however the dominant residue in this position is a glutamic acid. The 3D view of the complex of the neurogenic differentiation factor 1 (mouse homolog of TWIST1) and the transcription factor E2-alpha bound to DNA (PDB code 2QL2). The conserved and the specificity residues derived from multiple sequence alignment of TWIST1 homologs are mapped on the 3D structure of the neurogenic differentiation factor 1 and shown in blue and orange, respectively; DNA fragment is shown in violet; an alphahelix fragment of E2-alpha is shown in green. The glutamic acid E147 in the position of the mutated residue (shown in dark red) contributes to the complex formation by interacting with the charged residues of the arginine 600 and the main chain of the glutamic acid 601 of E2-alpha protein; it also interacts with the hydroxyl group of serine 138 of TWIST1. The mutation of the glutamine Q147 to glutamic acid in TWIST1 is predicted to affect the stability, the complex and the regulation of transcription performed by TWIST1 molecule.

Fig. 1 .

 1 Fig. 1. Sequencing chromatogram of novel variants identified in this study. The arrows indicate the positions of the variants. (Left) TWIST1 c.439C>G; (Right) FGFR1 c.730_731insG.

Fig. 2 .

 2 Fig. 2.

Fig. 3 .

 3 Fig. 3.

Table 1

 1 Syndromic Craniosynostosis with Sagittal Suture Involvement and Related Genes

	Reference	21		22		23	24		25		26			27	28	10,29	30	31	32		33		34		35			36	37	38	39
	Mechanisms/Pathway	Reduced ATP-dependent K+ channel	inhibition	Abnormal ALP enzyme expression	and mineralization	Affect GAG sulfation	Abnormal kinase activity and	secretion in mineralization	TGFβ		FGF			FGF	FGF		FGF	Disrupt GLI-PTCH1-SHH signaling	Reduced cilia frequency and length		Cilia		Misexpression and/or overexpression	of IHH signaling	Abnormal IL11/gp130 affect JAK/	STAT1/3 or SHP2/MAPK/ERK	phosphorylation cascades	Loss of Notch signaling activity	Perturb HH and Nodal signaling	FGF/BMP	Abnormal osteoclast function	Aberrant function in osteoblast-	mediated bone mineralization
	Mutation effect	LOF		LOF		LOF	LOF		LOF		11q11-q13.3 duplication			GOF	GOF		GOF	LOF	LOF		LOF		Microduplication at the	IHH locus	LOF			LOF	LOF	GOF	LOF	LOF
	Sutures affected	Coronal, sagittal		Coronal, sagittal, lambdoid		Sagittal	Multiple		Coronal, sagittal, metopic		Sagittal, metopic,	lambdoid, temporal,	squamosal	Multiple	Sagittal	Sagittal	Coronal, sagittal, lambdoid	Coronal, sagittal, metopic	Sagittal		Sagittal		Coronal, sagittal, metopic,	lambdoid	Coronal, sagittal, metopic,	lambdoid		Coronal, sagittal, metopic,	lambdoid	Coronal, sagittal, metopic	Sagittal, metopic	Coronal, sagittal	Coronal, sagittal
	OMIM Syndrome/Phenotype	601439 Cantu syndrome		171760 Hypophosphatasia infantile type		603799 Autosomal recessive Larsen syndrome	611061 Raine syndrome (non-lethal osteosclerotic bone	dysplasia)	134797 Shprintzen-Goldberg syndrome (Marfanoid	craniosynostosis)	164950, 164980 Syndromic multiple craniosynostosis			136350 Osteoglophonic dysplasia	176943 Crouzonoid syndrome; Familial scaphocephaly	syndrome	134934 Thanatophoric dysplasia type II	165240 Greig cephalopolysyndactyly syndrome	606045 Cranioectodermal dysplasia 1 (Sensenbrenner	syndrome)	614068 Cranioectodermal dysplasia 3 (Sensenbrenner	syndrome)	600726 Craniosynostosis and syndactyly		600939 Craniosynostosis and dental anomalies syndrome			601920 Alagille syndrome	604267 Carpenter syndrome 2	123101 Boston-type craniosynostosis	607649 Infant osteopetrosis craniosynostosis, Chiari	malformation type I	300550 X-linked dominant hypophosphatemic rickets
	Gene	ABCC9		ALPL		CHST3	FAM20C		FBN1		FGF3, FGF4			FGFR1	FGFR2		FGFR3	GLI3	IFT122		IFT43		IHH		IL11RA			JAGGED1	MEGF8	MSX2	OSTM1	PHEX

Table 2

 2 Mutations/Variants Found in Human Sagittal Nonsyndromic Craniosynostosis

	Reference	48		49	16			This study	15	50			11	10	20	51			52		53		16			14				This study
	Population Frequency (1KG/ESP/ExAC)	1/12956 (ESP); 1/49772 (ExAC)	0	0	21/5008 (1KG); 1/6788 (ESP); 85/35096	(ExAC)		0	0	0			1/121146 (ExAC)	0	1/13006 (ESP)	2/12988 (ESP); 4/121408 (ExAC)	24/5008 (1KG); 1/12988 (ESP); 247/121406	(ExAC)	0	2/121380 (ExAC)	1/5008 (1KG); 8/12988 (ESP); 129/121130	(ExAC)	0			0				1/120318 (ExAC)
	Mutation Detection Rate	1% (2/203)		1% (1/70)	1% (1/83)			1% (1/93)	2 patients	1 patient			3% (1/29)	1% (1/94)	3 patients	2% (2/94)			0.5% (2/431)		1% (1/104)		1% (1/83)			2 patients from 1 family, 1	has scaphocephaly, 1 has no	synostosis (reduced	penetrance)	1% (1/93)
	Nucleotide Change Amino Acid Change Case Population Studied	c.19G>T V7F 203 NSC patients (197 sNSC and 6	c.631A>G K211E multiple suture CS with sagittal suture involvement)	c.194G>A R65C 70 sNSC	c.1508C>T T261M 164 infants with isolated, single-suture	craniosynostosis (83 sagittal, 37 metopic,	12 lambdoid, 32 unilateral coronal)	c.730_731insG A244fs*26 93 sNSC patients	c.1032G>A A344A 2 familial sNSC	c.833G>A C278Y 1 familial (1 individual has isolated	sagittal/unilambdoid synostosis, 3 others	unaffected)	c.943G>A A315T 29 sNSC patients	c.1139_1196dup M400fs*53 630 craniosynostosis (94 sNSC)	c.1000G>A A334T 1 familial, sagittal and metopic ridging	c.1217G>A R406H 186 NSC (94 sagittal, 46 coronal, 46	c.2570A>G N857S metopic)		c.1346G>C S494T 431 NSC patients (no information if only	c.1640G>A C592Y sagittal or multiple sutures affected)	c.263G>A D68N 104 sagittal/coronal (no information if	nonsyndromic or multiple sutures)	c.563C>T S188L 164 infants with isolated, single-suture	craniosynostosis (83 sagittal, 37 metopic,	12 lambdoid, 32 unilateral coronal)	c.602C>A S201Y 178 patients (84 Saethre-Chotzen	syndrome, 7 Muenke syndrome, 12 other	syndromes, 75 nonsyndromic coronal	synostosis)	c.439C>G Q147E 93 sNSC patients
	Gene	AlX4		ERF	FGFR1				FGFR2						FGFR3	IGF1R			LRIT3		SOX6		TWIST1

Plast Reconstr Surg. Author manuscript; available in PMC 2017 March 01.

Plast Reconstr Surg. Author manuscript; available in PMC 2017 March 01.

Acknowledgments

We thank the patients and their families for participation in this study. We are indebted to Charlotte M. Druschel, M.D., M.P.H., who is retired from the Congenital Malformation Registry, New York State Department of Health, Albany, New York and the Department of Epidemiology and Biostatistics, School of Public Health, University at Albany State University of New York, Albany, New York. She was a major contributor to this project in design and recruitment of New York State participants. This research was supported by grants including CDC 5 R01 DD000350 and NIH/NICHD 1 P01 HD078233.

designed the study, coordinated and oversaw the project.

Statement of Ethical approval

Written informed consents were obtained from the parents or guardians of the patient. The study protocol was approved by the Icahn School of Medicine at Mount Sinai and the University of Iowa and New York State Department of Health institutional review boards, and conformed to the guidelines set forth by the Declaration of Helsinki.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.


Filipe M M Marreiros
email: filipe.marreiros@liu.se

Yann Heuzé
email: yann.heuze@u-bordeaux.fr

Michael Verius
email: michael.verius@i-med.ac.at

Claudia Unterhofer
email: claudia.unterhofer@tirol-kliniken.at

Wolfgang Freysinger
email: wolfgang.freysinger@i-med.ac.at

Wolfgang Recheis
email: wolfgang.recheis@i-med.ac.at

Custom implant design for large cranial defects

Keywords: Cranial reconstruction, reconstructive surgery, geometric morphometrics, radial basis functions and thin plate spline

The aim of this work was to introduce a Computer Aided Design (CAD) tool that enables the design of large skull defects (> 100cm 2 ) implants. Functional and aesthetically correct custom implants are extremely important for

patients with large cranial defects. For these cases pre-operative fabrication

Introduction

Functional and aesthetically correct custom implants are extremely important for patients with large cranial defects (> 100cm 2 ). Usually, small defects are treated using autograft material due to the high rate of a successful incorporation by defect site tissues. However, problems of donor site morbidity, sufficient donor material, bone tissue quality and crafting the correct shape from the available donor increase dramatically with the size of the defect. In these cases the pre-operative fabrication of an implant should be considered.

Using current medical imaging modalities, in our case 3D Computed Tomography (CT) and specialized software tools the anatomically correct design of such implants is possible [START_REF] Parthasarathy | 3D modeling, custom implants and its future perspectives in craniofacial surgery[END_REF], [START_REF] Liao | Using threedimensional multigrid-based snake and multiresolution image registration for reconstruction of cranial defect[END_REF], [START_REF] Rotaru | Cranioplasty With Custom-Made Implants: Analyzing the Cases of 10 Patients[END_REF], [START_REF] Liao | Three-dimensional reconstruction of cranial defect using active contour model and image registration[END_REF], [START_REF] Lüthi | Building shape models from lousy data[END_REF], [START_REF] Wu | Anatomically constrained deformation for design of cranial implant: Methodology and validation[END_REF], [START_REF] Hierl | CAD-CAM implants in esthetic and reconstructive craniofacial surgery[END_REF], [START_REF] Dean | Computer aided design of pre-fabricated cranial plates[END_REF], [START_REF] Min | Highly accurate CAD tools for cranial implants[END_REF], [START_REF] Min | Computer Aided Design of Cranial Implants Using Deformable Templates[END_REF], [START_REF] Eufinger | Computer-assisted prefabrication of individual craniofacial implants[END_REF], [START_REF] Hsu | Application of three-dimensional orthogonal neural network to craniomaxillary reconstruction[END_REF], [START_REF] Hsu | Application of three-dimensional orthogonal neural network to craniomaxillary reconstruction[END_REF] [START_REF] Carr | Surface interpolation with radial basis functions for medical imaging[END_REF], [START_REF] Wehmöller | CAD by processing of computed tomography data and CAM of individually designed prostheses[END_REF]. The different software tools vary in the user interaction required to design the implants: manual, semi-automatic or automatic. Semi-automatic or automatic methods aid considerably the design of such implants and a complete and coherent software CAD (Computer Aided Design) pipeline is of great importance to the user in charge of the production of the implants.

In the next sections we will present the pipeline blocks and methods of our CAD implant design tool that allows the generation of 3D virtual implants that can be further used to produce the physical implants.

Materials and methods

Several CAD tools exist that can be used for implant design ranging from manual to automatic solutions. The core of our CAD tool is a semi-automatic method that relies on the interpolation properties of the Radial Basis Function (RBF), in particular the Thin Plate Spline (TPS), and the statistical tools of Geometric Morphometrics (GM) [START_REF] Bookstein | Morphometric Tools for Landmark Data: Geometry and Biology[END_REF], [START_REF] Slice | Modern Morphometrics in Physical Anthropology[END_REF]. The GM methods belong to the more general statistical shape analysis methods, which is the analysis of the geometrical properties of sets of shapes by statistical methods. An introduction to these methods is provided in [START_REF] Stegmann | A Brief Introduction to Statistical Shape Analysis[END_REF]. Before we describe in detail the blocks that compose the CAD tool pipeline we will introduce the most important methods and compare to related work. All the defects presented in this work, including simulated defects used for evaluation are large skull defects (> 100cm 2 ).

Radial Basis Function

RBFs can be used as surface interpolators, the reconstruction properties of the RBFs have been described in detail by Carr et al. [START_REF] Carr | Surface interpolation with radial basis functions for medical imaging[END_REF], [START_REF] Carr | Reconstruction and representation of 3D objects with radial basis functions[END_REF]; the present work is inspired by these methods, and builds upon our previous work [START_REF] Heuzé | The use of procrustes average shape in the design of custom implant surface for large skull defects[END_REF]. RBFs are used in GM to study shape intersubject differences using deformation grids. In Computer Graphics (CG) RBFs are mainly used for interpolation purposes to generate smooth functions passing through input data (centers) and can be used, for instance, for hole filling.

In GM however, RBFs are used to encode differences between shapes, or to map one shape onto another. TPS warping and spline relaxation are among the most used applications of RBFs in GM.

The highly complex structure of the human skull in the facial area exhibits sharp edges and holes so that the presented method can only be used for the neurocranium, where TPS can be employed to minimize bending energy [START_REF] Bookstein | Principal warps: Thin-plate splines and the decomposition of deformations[END_REF].

GM -Landmark

A detailed landmark classification can be found in the work of Bookstein [START_REF] Bookstein | Morphometric Tools for Landmark Data: Geometry and Biology[END_REF]. In brief, semilandmarks are points on a structure allowing the incorporation of information from ridge curves and surfaces of the structure, i.e. a skull. Classical anatomical landmarks bear full threedimensional information, whereas semilandmarks present "deficient" coordinates and are computed on base of the available coordinates. Semilandmarks are treated as missing data and are estimated in order to minimize the bending energy between a template shape and the target shape. 2D semilandmarks and spline relaxation along curves were introduced by Bookstein [START_REF] Bookstein | Landmark methods for forms without landmarks: morphometrics of group differences in outline shape[END_REF], work that was extended to three dimensions by Gunz et al. [START_REF] Gunz | Semilandmarks in three dimensions[END_REF]. A curve, surface or fully relaxed semilandmark is a landmark with one, two or three unknown dimensions and thus semilandmarks can slide on curves, surfaces and in 3D space. Anatomical landmarks are not relaxed, they are fully defined.

Fully relaxed semilandmarks are generally used to estimate semilandmarks of regions where missing data occur, i. e. defect sites. We will not differentiate between landmarks and semilandmarks, both will be distinguished by relaxation type: anatomical, curve relaxed, surface relaxed and fully relaxed landmark. Missing landmarks are any of the previous ones that are used to delineate the defect region for implant placement with a slight overlap to correctly compute the implant boundary. In most cases the complete set of missing landmarks is formed by surface relaxed and fully relaxed landmarks.

GM -methods used

From the suite of GM methods we use the Procrustes Average Shape (PAS), TPS spline relaxation and Generalized Procrustes Analysis (GPA) [START_REF] Rohlf | Extensions of the procrustes method for the optimal superimposition of landmarks[END_REF]. The PAS is computed using a database of 130 anonymized CT exams with axial slices and average slice thickness of 1 mm. Demographically, 73 males (average age 51.14 years, standard deviation 16.28 years, maximum 82 years, minimum 21 years) and 57 females (average age 57.26 years, standard deviation 21.32 years, maximum 97 years, minimum 20 years) are in the data set. PAS and GPA align the population of shapes iteratively and are well known [START_REF] Gunz | Semilandmarks in three dimensions[END_REF], [START_REF] Heuzé | The use of procrustes average shape in the design of custom implant surface for large skull defects[END_REF]. In every iteration each specimen is relaxed to the PAS and is aligned using GPA. With all shapes aligned and relaxed we compute the PAS for the actual iteration. The final PAS should be taken from the final iteration.

The results of this operation are: the PAS (that might happen to be the best fitting shape) and all landmark sets in the database relaxed to the PAS. To compute the best fitting shape we also make use of GPA, for more details see sections: measuring the anatomical landmarks, initial landmark estimation and fully relaxed landmark selection. The fully relaxed landmarks are obtained using symmetry information (if possible) and the best fitting shape.

Related work

The closest work to ours is by Min and Dean [START_REF] Dean | Computer aided design of pre-fabricated cranial plates[END_REF], [START_REF] Min | Highly accurate CAD tools for cranial implants[END_REF], [START_REF] Min | Computer Aided Design of Cranial Implants Using Deformable Templates[END_REF], who also use landmarks and a TPS warp [START_REF] Bookstein | Landmark methods for forms without landmarks: morphometrics of group differences in outline shape[END_REF] for the reconstruction. Despite a basically similar concept, our approach differs particularly with respect to the determination of the defect boundary, implant triangle mesh generation and the number of specimens that provide the estimate for missing landmarks. The present approach, however, goes beyond their work by incorporating symmetry and the best fitting shape instead of either symmetry or mean shape. We also used the TPS spline relaxation instead of the TPS warp providing a more accurate positioning of the landmarks, as explained later on in the text. The last main difference is that we developed a method to find the best fitting shape from our database; in this paper we will show that using the best fitting shape combined with symmetry is in most cases superior to just using a mean shape.

The main difference to Carr's surface interpolation with RBFs for medical imaging [START_REF] Carr | Surface interpolation with radial basis functions for medical imaging[END_REF], is that the missing landmarks and centers of the RBF, of the defect site are obtained by GM methods. Carr's surface interpolation uses the vertices of the defect's boundary region as the centers of the RBFs.

Wu et al. [START_REF] Wu | Anatomically constrained deformation for design of cranial implant: Methodology and validation[END_REF] used anatomical constrained deformations based on reference models (3D-matching and adaptive deformation) and compared their method against RBFs. They found that RBFs provided worse results as only the boundary information was used as centers for the RBFs.

Finaly, Lüthi et al. [START_REF] Lüthi | Building shape models from lousy data[END_REF] showed that Procrustes alignment is used to generate a statistical shape model using probabilistic principal component analysis and Liao et al. [START_REF] Liao | Using threedimensional multigrid-based snake and multiresolution image registration for reconstruction of cranial defect[END_REF], [START_REF] Liao | Three-dimensional reconstruction of cranial defect using active contour model and image registration[END_REF] used active contours and image registration for the reconstruction. These methods can be used also to estimate missing data; however, they do not allow to incorporate symmetry information, one of the key sources of information of our method.

CAD design pipeline

The CAD tool is built upon the Visualization Toolkit (VTK) [START_REF] Schroeder | The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics[END_REF] and the design pipeline is shown in Figure 1. The process starts by performing a manual threshold segmentation (a value of 250 was used in all cases) of the 3D CT DICOM data. Since the region of interest is the largest connected bony region, all small bone fragments and artifacts have to be removed. The next step is to manually measure the anatomical landmarks. Followed by an initial estimate of all remaining landmark positions. In the defect region the fully relaxed landmarks are estimated using the TPS relaxation. After n relaxation iterations the missing landmarks are manually selected and used for generating the outer surface mesh of the implant. The outer mesh serves to calculate the mesh of the implant's inner surface. Furthermore the surface has to be fitted to the defect boundary. This implies identifying the triangles of the defect boundary and finally to connect both the inner and the outer meshes to the final implant mesh.

Removal of artifacts and fragments

Artifacts and fragments removal is done with the vtkImageSeedConnectivity, which is a connected component labeling algorithm [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF], [START_REF] Shapiro | Computer Vision[END_REF]. Since we know that the CTs contain only skulls, and that the largest connected region in the (multiplanar reformatted) medial sagittal plane corresponds to the skull and not a fragment or an artifact, we can on this slice test every voxel above threshold for its connectivity check the largest connectivity region. The pixels (actually voxels) with the largest connected region are selected for a 3D volume connectivity test. Finally, the largest connected region is kept and all other removed.

Measuring the anatomical landmarks, initial landmark estimation and fully relaxed landmark selection

17 cranial anatomical landmarks are measured; some are standard landmarks for GM studies [START_REF] Ousely | Data Collection Procedures for Forensic Skeletal Material[END_REF], [START_REF] Aeillo | An introduction to human evolutionary anatomy[END_REF], all are listed in Table 1. Two sets of standard landmarks are non-standard: [START_REF] Dean | Computer aided design of pre-fabricated cranial plates[END_REF][START_REF] Min | Highly accurate CAD tools for cranial implants[END_REF] and [START_REF] Hsu | Application of three-dimensional orthogonal neural network to craniomaxillary reconstruction[END_REF][START_REF] Hsu | Application of three-dimensional orthogonal neural network to craniomaxillary reconstruction[END_REF], were introduced to form a more complete landmark distribution around the skull. A ray-casting approach [START_REF] Drebin | Volume rendering[END_REF], [START_REF] Levoy | Display of surfaces from volume data[END_REF], [START_REF] Levoy | Efficient ray tracing of volume data[END_REF] was used to measure/set the landmarks via a trilinear interpolated, ray sampled at 0.05 mm. Figure 2 shows some of the anatomical landmarks on the skull. An initial estimate of the semilandmarks was obtained applying different projection schemes, by using different sets of anatomical landmarks (according to the anatomical location) to generate projection planes or lines where rays are casted to the surface. This is exemplified on one of these regions, the region defined by the bilateral landmarks Fronto Temporale and Zygomatic Root. These landmarks define a plane that allows defining circular sampling lines perpendicular to it with the respective center at midline, defined by the line connecting the centroids of the two pairs of bilateral landmarks. For demonstration purposes a few sampling and initialization points on the boundary (green and blue points, respectively) are shown. This ray-casting is defined via the projection centers along the midline (red points in Figure 3) and the corresponding circular sampling points (green dots in Figure 3). Corresponding (surface relaxed) landmarks are generated if the ray reaches the last voxel above or equal to a threshold; else the landmark is considered fully relaxed. The initial estimates for a complete skull (727 landmarks) are shown in Figure 4. Missing anatomical landmarks need to be estimated for defective skulls. All missing landmarks can be estimated via symmetry, or else via the best fitting shape (defined below). To estimate these anatomical landmarks via symmetry we first calculate the side (left or right) of the skull that has most fully relaxed landmarks, this will be defined as target and the opposite side as the template. The bilateral missing landmarks are excluded. The template set will contain all the defined landmarks of the corresponding side while the target set the symmetric points, if these are missing they are set to fully relaxed. In this case the landmarks do not slide (explained below -only one iteration is required). After the TPS relaxation the fully relaxed landmarks will hold the new landmark positions. This procedure is repeated twice, for each side of the skull. For the bilateral missing landmarks the best fitting shape is used, by setting the missing land-marks to fully relaxed and perform TPS relaxation to obtain their position. Once the anatomical landmarks are known the remaining landmarks can also be obtained in a similar fashion. Fig. 4 Landmarks placed on the skull surface using the initial estimation method (727 landmarks). The initial estimation uses anatomical landmarks (cyan) to define projection schemes used to place the surface landmarks (violet).

In the symmetry case the TPS is used considering that the available landmarks (surface landmarks) are not relaxed, they do not slide, because surface normals are pointing generally in opposite directions for bilateral landmarks, this means that the tangent planes are not equivalent consequently the sliding is not possible. Once we obtain the landmark estimation of one side we have to repeat the operation for the remaining side.

The best fitting shape is obtained via GPA alignment of the patient shape and all relaxed shapes in the database including the PAS without considering the fully relaxed landmarks. The set (patient shape, database shape or PAS) with the smallest Root Mean Square Error (RMSE) between the landmark positions of the patient shape and each shape in the database is taken as the best fitting shape. Some landmarks still exist that have to be set manually, c. f. Figure 5, e. g. at the inner surface of the skull or at the boundary region of the defect.

Fig. 5 Initial landmark estimation of a defective skull. Landmarks on the inner surface of the skull and on the boundary region of the defect need to be selected manually to full relaxation (arrows point to some of these landmarks).

TPS relaxation

The TPS relaxation is used to obtain the PAS, to estimate missing data and to slide surface landmarks to positions with minimum bending energy. All skulls in the database are complete (with no defect) with anatomical and surface landmarks. For the relaxation calculations the anatomical landmarks are considered as surface landmarks. Setting anatomical landmarks to surface relaxed landmarks will allow sliding all the landmarks to positions providing minimal bending energy between target and template shapes. Moreover, user subjectivity in landmark measurement is eliminated. Some TPS generated landmarks might be slightly misplaced with respect to the correct anatomical sites (mainly due to the user initial setting); however, they are guaranteed to be in homologous positions.

For analyzing defective skull data the anatomical landmarks are also set to surface relaxed landmarks. All available surface landmarks are used and the best fitting shape is iteratively searched.

GM is frequently performed with TPS relaxation on triangular meshes, but never, as far as we know, was the voxel information used directly. The process is started by computing the surface relaxation tangent plane were the landmarks can slide. The normal of the plane can be computed with a 3D extension of the 2D Sobel operator [START_REF] Sobel | A 3 x 3 isotropic gradient operator for image processing[END_REF] that uses trilinear interpolation to generate isotropic data. The plane can be defined using the normal and the landmark position. Some tangent planes and the normal vectors at the landmarks of a complete skull are displayed in Figure 6. Landmark relaxation slides the landmarks on the tangent planes so that they eventually move off the anatomic surface under consideration and need to be projected back on the surface. It is essential for the whole process that this projection generates landmarks on the outer surface. This is assured by sampling around the landmarks in a 101x101x101 voxels bounding box to classify whether voxels are in the outer surface or not: its value has to be at least the bone threshold and at least one of its 26 neighbors has to be below threshold. All other possibilities are classified as non-boundary voxels, ideally providing two sets of voxels corresponding to inner and outer surfaces, respectively. For the following, the outer surface only is considered. The new landmark position will correspond to the voxel with the smallest distance to the previous landmark position.

Exterior surface interpolation using 2D RBF

The interpolation of the exterior surface is done with a 2D RBF [START_REF] Carr | Surface interpolation with radial basis functions for medical imaging[END_REF] that allows the construction of a depth-map centered at the (manually selected) missing landmarks of the defective skull. In order to close the implant boundary to the rest of the rest of the defective skull, some of the missing landmarks are placed on the bone, see later.

The depth-map requires a plane that allows the correct implant reconstruction. This plane is defined by the normal vector chosen as the actual direction of the surgical implant insertion process and a point at the centroid of the missing landmarks. One point of this vector is defined as the centroid of all the landmarks, the other is the centroid of the missing landmarks. The exterior surface (Figure 7) is obtained by evaluating the depth-map at points corresponding to the vertices of the mesh on a regular triangle grid (z = 0) with adaptable numbers of columns and rows in order to create triangles greater than the largest neighbor voxel distance.

Triangle classification and connectivity test

The connectivity test is crucial for implant design. The mesh is sampled regularly along the triangle surface and interpolated trilinearly using its eight neighbor voxels. Each triangle is labeled: "1" In space: All sampling gray values are below the threshold. "2" Inside the bone: All sampling gray values are above the threshold. "3" Else: Mixed. All triangles with the label "2" intercept the outer surface of the skull and can be discarded. The implant corresponds to the largest connected region formed by triangles with label 1(cyan region, Figure 8) that can be chosen simply by counting its triangles on the regular grid.

In figure (Figure 8) the two possible test configurations [START_REF] Deutsch | Thinning algorithms on rectangular, hexagonal, and triangular arrays[END_REF] with a maximum of twelve connected neighbors (gray triangles) are presented, the configuration to be used is dependent on the triangle type: upper or lower triangle.

This test, similar to the connected components labeling algorithm [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF], [START_REF] Shapiro | Computer Vision[END_REF], is applied to the triangles, if no neighbor triangles share the value, a new group is formed. If only one shares the label, it is added to this labeled group; if more Fig. 7 Exterior mesh obtained using the 2D RBF. The yellow points are missing landmarks, the other points surface landmarks, the surface is shown in cyan.

than one connected triangles share the same label a new group is created and groups are merged.

Boundary fitting

The boundary triangles (red triangles, Figure 8) that have value "3" are directly connected to the triangles of the largest connected region. Triangles sharing only one edge with (i) the remaining boundary triangles and (ii) the triangles of the largest connected region are removed and now all boundary triangles have at least two connecting edges, and the fitting process can be initialized.

In order to fit the implant to the defect region, it is necessary to ensure that its boundary triangles are completely outside the bone. The process starts by independently sampling along the inner edges of border triangles of the implant as depicted in Figure 10 (top). The sampling is stopped when a sampling position below threshold (outside of bone) is found; at the previous step there is still some overlap (Figure 10, top right). Next, the border edge is sampled Figure 10 (bottom, left), if the edge has "inside sampling points" then the edge get discarded and the triangle vertices of the new edge are updated, this process is repeated until one edge with only "outside sampling points" is found. Furthermore, in order to produce a small gap between bone and Fig. 8 Outer surface mesh as the largest connected region (selection of triangles with value "1". The red triangles are directly connected to the cyan region and have value "3" (boundary triangles). The yellow points and violet points are fully relaxed and surface relaxed landmarks, respectively. (Bottom) Sampling the non connected edge (several iterations might be performed, updating the edge position, until an edge containing only outside sampling points is found). The gray area represents the bone defect area, the red dots represent the sampling points inside, and the blue dots represent the points outside. The orange arrows represent the direction of the edges sampling and the small green arrows represent the sampling performed with a direction parallel to the implant vector, this is used to provide a small gap between bone and implant, in our case a maximum distance of 1 mm from the origin sampling point was used. If the sampling positions along these directions contain inside sampling points then the correspondent origin sampling point is considered also an inside sampling point, otherwise an outside sampling point.

implant, sampling was performed with directions parallel to the implant vector (Figure 10 -green arrows), in our case a maximum distance of 1 mm from the origin sampling point was used, but this can be chosen freely to control the gap size that allows accounting for accumulated errors and ease intraoperative placement. If the sampling positions along these directions contain "inside sampling points" then the correspondent origin sampling point is considered also an "inside sampling point", otherwise an "outside sampling point".

Creating the implant interior surface and its connection onto the exterior surface

To correctly represent the inner profile of the boundary region a much higher triangle density than on the outer surface is needed.

Figure 12 shows the two-dimensional sketch of the effects of different sampling rates. The inner surface of the implant does not resemble the outer surface of the bone at locations where the implant should fit: red lines (implant) are seen to intersect the bone. This can be attributed to insufficient sampling [START_REF] Nyquist | Certain topics in telegraph transmission theory[END_REF], [START_REF] Shannon | Communication in the presence of noise[END_REF]. We use a triangle subdivision scheme [START_REF] Jin | Subdivision interpolating implicit surfaces[END_REF] for the outer surface that stops if a maximum set of divisions or if all triangles edges are shorter than the minimum voxel edge distance. The implant outer surface vertices are ray-casted towards the bone surface along rays parallel to the implant vector direction Figure 12, similar to [START_REF] Drebin | Volume rendering[END_REF], [START_REF] Levoy | Display of surfaces from volume data[END_REF], [START_REF] Levoy | Efficient ray tracing of volume data[END_REF]. A maximum implant thickness can be set. For safety purposes the rays are pushed back by 0.2 mm from their stopping points.

The last step is to connect inner and outer surfaces. First the outer boundary edge is broken up to match the number of vertices on the inner surface with the help of the previous process and the vertices are connected correctly (Figure 13).

Evaluation

The TPS spline relaxation is known to converge delivering sets of landmarks that minimize bending energy between the shapes involved. We used the evolution of the RMSE between iterations as an indicator of convergence. This test was applied to all skulls in the database, obtaining a combined RMSE of all possible landmark sets. And also the RMSE values for one complete skull wich is not from the database, one skull with an unilateral defect (but with all anatomical landmarks) and one skull with an bilateral defect (three anatomical landmarks are missing). The PAS was not tested because it depends on the alignment of all shapes in the database producing a slightly different PAS pose at every iteration. Since all shapes of the database are being tested the PAS test would be redundant, because if all shapes converge also the mean shape should converge.

The next test was an evaluation of the landmark position after relaxation, using four different approaches, namely testing symmetry, best fitting shape and PAS. The combinations are: Five complete skulls, not contained in the database were used. Surface landmarks were manually selected and converted to fully relaxed landmarks at different regions of the skull to simulate defects. Ideally these landmarks should be placed on the skull surface after relaxation. We computed the mean error (µ), maximum error (max) and percentage of landmarks with error less or equal to 1 mm. The last metric was chosen because the most similar previous work have used it, enabling a comparison. The errors were calculated for fully relaxed landmarks only, and were measured as the distance from the estimated positions to the closest point on the exterior skull surface (ground truth). Ten iterations of the TPS relaxation were performed. Region A includes the parietal and occipital bones; region B the frontal, sphenoid and temporal bones; region C the frontal and parietal bones; region D the frontal bone(including one orbit) and finally region E the frontal bone (including both orbits). An example of one artificial skull defect (yellow landmarks are fully relaxed) is shown in Figure 14 corresponding to Specimen 2 region A.

Results

The results of the convergence test can be seen in Table 2. In all cases after approximately 10 iterations the RMSE drops to considerable small values, thus the method is converging correctly, with the calculated landmark sets properly representing the shape of the skull region under consideration, and 10 iterations seem to be sufficient. The results of the landmark position evaluation, after relaxation are presented in Figure 15. The overall means are presented in Table 3. From the results we can observe that when using symmetry information a higher accuracy was achieved, but there was no significant difference between using the Best Fitting Shape or the PAS. The values of the four tests presented, prove the robustness of this method, even in the worst case (no symmetry only PAS -Region E, Specimen 3) a mean error of 1.01 mm, maximum error 4.39 mm, and 66.17% of the fully relaxed landmarks have an error less or equal to 1 mm. From all the cases studied the maximum mean error was approximately 1 mm, a value satisfactory for the clinical problem we are considering.

Discussion

The TPS methods require landmarks that are initially placed in anatomical locations. TPS relaxation is used to regularize the landmark positions by sliding them on the surface of the skull and using the bending energy to guide the relaxation. After the relaxation the landmark positions can be compared to the landmarks of the skulls in the database. As shown in Table 2 the landmark sets within a few relaxation iterations have minimum variation between iterations and the calculated landmark sets properly represent the shape of the skull region under consideration. The small differences that are still present in the last iterations can mainly be attributed to CT resolution affecting the surface tangents.

The values of the four tests presented in Figure 15, prove the robustness of this method, even in the worst case (no symmetry only PAS -Region E (frontal bone), Specimen 3) a mean error of 1.01 mm, maximum error 4.39 mm, and 66.17% of the fully relaxed landmarks with an error less or equal to 1 mm is found. Some regions provide better estimations, region A (parietal and occipital bones) and region E (frontal bone) are worst. The reason region A produces less accurate results is due to a smaller density of landmarks in comparison to other regions; this should be corrected in the future. Regarding region E it is expected to be the most difficult region to estimate [START_REF] Poukens | A classification of cranial implants based on the degree of difficulty in computer design and manufacture[END_REF].

It can be expected that symmetry and the best fitting shape provide the best solutions. But as Table 3 shows the values obtained are almost similar to using symmetry and PAS. This might be due to a suboptimal population of our database. Our results proof that symmetry provides better estimates. The complete skulls used for the evaluation process were not particularly asymmetrical, they are representative for the normal population, though exceptions could of course exist especially when the patient suffered trauma or pathologies producing some levels of asymmetry. Our method is able to compensate to some extent for asymmetric cases because it incorporates information of the landmarks surrounding the defect, this way constraining the process to obtain the missing landmarks. Furthermore, for the symmetry cases there are always two sets of points Source and Target (from opposite sides of the skull -each landmark is labeled either left of right). The missing landmarks are estimated using the TPS that uses the least bending energy to deform these shapes, i.e. a combination of information from both sides is used. This is in fact the main difference between our method and traditional symmetry methods that typically only use information from the opposite side of the defect and try to fit it to the defect side.

Wu et. al [START_REF] Wu | Anatomically constrained deformation for design of cranial implant: Methodology and validation[END_REF] reported an average of 81.6% (68% -96%) surface points with an error less than 1 voxel. In our case using symmetry and the best fitting shape an average of 89.29% (72.64% -100%) missing landmarks with an error less or equal to 1 mm was obtained. Wu's method used a mixture of small and large defects and reported for the large defects an average of 80.5% surface points within 1 voxel error.

When compared to symmetry methods like (mirrored volume) our method can look more complicated. From a user perspective with our method the user is required to identify anatomical landmarks (17 points) and once the landmark initial projections have been performed it is also necessary to select the missing landmarks; while in the mirror volume techniques the user must try to fit the mirrored volume to the defect area and most likely perform several manual post processing steps before the implant shape is finally obtained. It is crucial here to highlight the user-dependent dimension of such manual post processing and the related poor reproducibility and repeatability of such method, not to mention the time required to perform such adjustments. Though this manual approach might look less complex, we think that this approach can reveal itself as complex and time consuming without providing systematic robust results. Ultimately, if the defect is bilateral the mirror volume cannot be used.

The computation of the implant vector could be performed using alternative methods for irregularly distributed landmarks. But the present approach is simple and suffices for the majority of cases and allows for manual correction.

The boundary fitting method is not intended to produce the most accurate approximation. The goal, instead, is to obtain a smooth boundary free from sharp edges that follows the main contours of the boundary with a surgical safe margin.

All applied methods were designed for the 2D RBF case, including the connectivity optimizations. One might argue to use 3D RBFs instead [START_REF] Carr | Reconstruction and representation of 3D objects with radial basis functions[END_REF], [START_REF] Turk | Modelling with implicit sufaces that interpolate[END_REF] but for the majority of cases investigated 2D RBFs are sufficient. Only for cases where ambiguities in the depth-map occur, most likely to happen at orbital regions, 3D RBFs will be the method of choice. Another advantage of the 3D RBFs is that missing landmarks do not have to be selected as the 3D RBFs provide a signed distance measure for testing whether each point (voxel) is inside, outside or on the boundary of the volume. With this signed distance measure all points can be classified. At this time, this method has two limitations. First, if the defect area is not or is insufficiently covered by the landmark set, the reconstruction process will not be ideal and may possibly fail. Second, if a leaking effect occurs during the connectivity test due to empty regions. Figure 18 shows such a case (the frontal bone forming the top of the orbit is part of the defect area). This problem can possibly be overcome by using 3D RBFs.

The present approach, however, is simple and computationally cheap and can also be used to reconstruct other bony regions with smooth surfaces such as the femur. In an anthropological context the method could be valuable for reconstructing fragmentary bone remains.

Conclusion

An original and innovative method for constructing implants to cover large unior bilateral defects to the neurocranium is presented by combining symmetry and best fitting shape.

The whole design process can be used as a CAD software that allows producing model files of the implant for fabrication.

Fig. 1

 1 Fig. 1 Implant design pipeline.

Fig. 2

 2 Fig. 2 Anatomical Landmarks (cyan) as seen from different views. The tubular structure is a catheter for ventricular drainage and is shown as it has a Hounsfield units similar to bone.

Fig. 3

 3 Fig. 3 Projection scheme: The solid line represents the plane formed by the 4 anatomical bilateral landmarks Fronto Temporale and Zigomatic Root (light blue). The dark blue points are sampling positions along the lateral boundaries used to define the semicircles orthogonal to the plane. The green points are sampling positions along the semicircles, the red points are sampled in the midline. The projection is started at the red points radially in the direction of the sampling points (arrows). The rays stop at the last voxel above or equal to threshold, providing the new landmark positions.

Fig. 6

 6 Fig. 6 Normals (green lines) of the local tangent planes (yellow) for each landmark.

Fig. 9

 9 Fig. 9 The two test configurations used for triangle connectivity, depending on the tupe of triangle tested: upper (bottom) or lower (top) triangle. The green circles are the corresponding vertices used to search for the connected neighbor triangles (gray).

Fig. 10

 10 Fig. 10 Placement of the boundary triangles outside the bone. (Top) Sampling along the connected edges.(Bottom) Sampling the non connected edge (several iterations might be performed, updating the edge position, until an edge containing only outside sampling points is found). The gray area represents the bone defect area, the red dots represent the sampling points inside, and the blue dots represent the points outside. The orange arrows represent the direction of the edges sampling and the small green arrows represent the sampling performed with a direction parallel to the implant vector, this is used to provide a small gap between bone and implant, in our case a maximum distance of 1 mm from the origin sampling point was used. If the sampling positions along these directions contain inside sampling points then the correspondent origin sampling point is considered also an inside sampling point, otherwise an outside sampling point.

Fig. 11

 11 Fig. 11 Exterior implant surface after boundary fitting.

Fig. 12

 12 Fig. 12 Projection of the outer surface at different sampling rates to obtain the inner surface.

Fig. 13

 13 Fig. 13 Connecting inner (green) and outer (blue) surfaces.

Fig. 14

 14 Fig. 14 Artificial skull defect (replacing surface landmarks (violet points) by fully relaxed landmarks (yellow points) and performing 10 TPS relaxation iterations). The region covered by the artificial defect involve the parietal and occipital bones.

1 .

 1 Symmetry and best fitting shape 2. Symmetry and PAS 3. No symmetry only best fitting shape 4. No symmetry only PAS

Fig. 15

 15 Fig.15Distance from estimated positions to the closest point on the exterior skull surface: Mean error (µ), maximum error (max), percentage of landmarks with error (≤ 1 mm). The calculations were performed with the results of the four tests: symmetry and best fitting shape, symmetry and PAS, no symmetry only best fitting shape and no symmetry only PAS. Ten iterations of the TPS relaxation were performed. Regions: A: parietal and occipital bones: B: frontal, sphenoid and temporal bones; C: frontal and parietal region; D: frontal region including one orbit; E: frontal region including both orbits. Five different specimens were tested. In the graphs (A1) corresponds to region A, specimen 1, the remaining cases follow the same convention.

Fig. 16

 16 Fig. 16 Bilateral defect implant.

Fig. 17

 17 Fig. 17Unilateral defect implant. Axially clipped to allow visualizing surface fitting.

  Fig. 17Unilateral defect implant. Axially clipped to allow visualizing surface fitting.

Fig. 18

 18 Fig. 18 Reconstruction limitation. A part of the frontal bone forming the top of the orbit is missing resulting in a leaking effect during the connectivity test.

Table 1

 1 Anatomical landmarks used in this study.

	IDs	Names
	1	Nasion
	2,3	Frontomalare Orbitale. (right, left)
	4,5	Frontomalare Temporale. (right, left)
	6,7	Fronto Temporale. (right, left)
	8,9	The most medial point in the sphenoid
		bone that is parallel to the frontomalare
		orbitale when the skull is viewed
		laterally. (right, left)
	10,11	Zygomatic Root (right, left)
	12,13	The point in the Zygomatic Process
		of the temporal bone most
		anterior medial. (right, left)
	14	Opisthion
	15	Basion
	16,17	Maxillo Frontale (right, left)

Table 2

 2 Landmark convergence

	Iteration	Combined RMSE	RMSE of complete	RMSE of skull	RMSE of skull
		of all skulls in	skull not in	with unilateral	with bilateral
		the database (mm)	the database (mm)	defect (mm)	defect (mm)
	1-0	3.08	3.04	-	-
	2-1	1.35	1.20	1.00	1.68
	3-2	0.80	0.83	0.57	1.04
	4-3	0.51	0.64	0.31	0.65
	5-4	0.37	0.35	0.22	0.44
	6-5	0.30	0.69	0.17	0.29
	7-6	0.26	0.59	0.13	0.26
	8-7	0.23	0.47	0.11	0.22
	9-8	0.22	0.39	0.11	0.23
	10-9	0.20	0.30	0.09	0.25
	11-10	0.19	0.18	0.10	0.24
	12-11	0.19	0.15	0.14	0.24
	13-12	0.19	0.12	0.10	0.21
	14-13	0.18	0.13	0.09	0.18
	15-14	0.17	0.11	0.08	0.20

Table 3

 3 Overall means of the: mean error (µ), maximum error (max) and percentage of landmarks with error ( ≤ 1mm) , using the four different methods.

		µ (mm) ≤ 1 (%) max (mm)
	Symmetry and Best Fitting Shape	0.42	89.29	1.79
	Symmetry and PAS	0.42	89.93	1.74
	No Symmetry only Best Fitting Shape	0.53	83.68	2.26
	No symmetry only PAS	0.52	85.05	2.24

Acknowledgements The authors would like to thank Brenda Frazier for editorial suggestions. We would like also to thank Philipp Gunz and Demetrios Halazonetis for their explanations of GM methods used in this work.

This work was supported by the EU FP6 Marie Curie Actions EVAN, contract number: MRTN-CT-2005-019564.

Conflict of interest

The authors declare that they have no conflict of interest.

Ethical approval

For this type of study formal consent is not required, since all the patient data was anonymized.