David Arthur Simpson
email: david.simpson@qub.ac.uk

Graeme Richard Clark

Sharon Alexander

Giuliana Silvestri

Colin Eric Willoughby

Molecular diagnosis for heterogeneous genetic diseases with targeted highthroughput DNA sequencing applied to retinitis pigmentosa

Keywords: Diagnostics test, clinical genetics, Genetic screening/counselling, Molecular genetics, Ophthalmology

Background

The genetic heterogeneity of many Mendelian disorders, such as retinitis pigmentosa which results from mutations in over 40 genes, is a major obstacle to obtaining a molecular diagnosis in clinical practice. Targeted high-throughput DNA sequencing offers a potential solution and was used to develop a molecular diagnostic screen for retinitis pigmentosa patients.

Methods

A custom sequence capture array was designed to target the coding regions of all known retinitis pigmentosa genes and used to enrich these sequences from DNA samples of five patients. Enriched DNA was subjected to high-throughput sequencing singly or in pools and sequence variants were identified by alignment of up to 10 million reads per sample to the normal reference sequence. Potential pathogenicity was assessed by functional predictions and frequency in controls.

Results and conclusions

Known homozygous PDE6B and compound heterozygous CRB1 mutations were detected in two patients. A novel homozygous missense mutation (c.2957A>T; p.N986I) in the cyclic nucleotide gated channel beta 1 (CNGB1) gene predicted to have a deleterious effect and absent in 720 control chromosomes was detected in one case in which conventional genetic screening had failed to detect mutations. The detection of known and novel retinitis pigmentosa mutations in this study establishes high-throughput DNA sequencing with DNA pooling as an effective diagnostic tool for heterogeneous genetic diseases.

INTRODUCTION

The completion of the Human Genome Project in 2003 was heralded as the dawn of an era of genomic medicine, [START_REF] Guttmacher | Welcome to the genomic era[END_REF] in which information from genomes would guide clinical decision making and deliver personalized medicine. [START_REF] Ginsburg | Genomic and personalized medicine: foundations and applications[END_REF] It was anticipated that accelerated detection of disease-related mutations would improve genetic diagnosis and prognosis. [START_REF] Burke | Genetic testing[END_REF] However, delivery of personalized genomic medicine requires not only access to the complete human genome, but availability of appropriate genetic tests for individual patients. The genetic heterogeneity of many Mendelian disorders is a major obstacle to obtaining molecular diagnoses in clinical practice. [START_REF] Tucker | Massively parallel sequencing: the next big thing in genetic medicine[END_REF] For example, Retinitis Pigmentosa [START_REF] Hartong | Retinitis Pigmentosa[END_REF] (RP[MIM #268000]), the commonest inherited retinal degeneration (prevalence 1:4000), is caused by mutations in over 40 genes for the non-syndromic form of the disease alone (Retnet: Retinal Information Network, http://www.sph.uth.tmc.edu/retnet/). Molecular genetic testing is important for clinical care, [START_REF] Drack | From the laboratory to the clinic: molecular genetic testing in pediatric ophthalmology[END_REF] enabling assignment of risk, genetic counseling and prognosis, and will be essential for enrolling patients in the future gene therapy trials likely to stem from the promising current human trials of RPE65 (MIM +180069) therapies for Leber congenital amaurosis. [START_REF] Cremers | Promises and challenges of genetic therapy for blindness[END_REF][START_REF] Smith | Prospects for retinal gene replacement therapy[END_REF] Apart from the genetic heterogeneity of RP, there are a number of obstacles that currently limit molecular diagnosis and therefore hinder the potential for personalized genomic medicine to guide clinical decision making. The lack of clearly defined genotype-phenotype correlations makes it difficult to direct testing to specific candidate genes. Approximately 50% of all RP patients are the only known affected family member and have no evidence of consanguinity; these patients are classified as isolated, sporadic or simplex cases [START_REF] Fishman | Retinitis pigmentosa. Genetic percentages[END_REF] which excludes diagnostic evaluation based on inheritance pattern. Traditional genetic screening for RP is laborious, although technological advances have had some impact. [START_REF] Koenekoop | Genetic testing for retinal dystrophies and dysfunctions: benefits, dilemmas and solutions[END_REF][START_REF] Mandal | Sequencing arrays for screening multiple genes associated with earlyonset human retinal degenerations on a high-throughput platform[END_REF][START_REF] Clark | Development of a Diagnostic Genetic Test for Simplex and Autosomal Recessive Retinitis Pigmentosa[END_REF] The development of massively parallel or 'next-generation' sequencing techniques which generate millions of DNA sequence reads in parallel during a single experimental run offers a potential solution. [START_REF] Tucker | Massively parallel sequencing: the next big thing in genetic medicine[END_REF][START_REF] Voelkerding | Next-generation sequencing: from basic research to diagnostics[END_REF] However, to date, applications of this high-throughput DNA sequencing approach as a molecular diagnostic tool have been limited due to the costs and perceived technical and data handling challenges. [START_REF] Vasta | Next generation sequence analysis for mitochondrial disorders[END_REF][START_REF] Chou | DNA sequence capture and enrichment by microarray followed by next-generation sequencing for targeted resequencing: neurofibromatosis type 1 gene as a model[END_REF] The aim of this work is to demonstrate that high-throughput DNA sequencing is now sufficiently established to warrant introduction into clinical practice. Herein we report the use of this technique to detect known and novel mutations in RP patients, demonstrating the clinical utility of this new technology.

METHODS

Participants

All applicable institutional and governmental regulations concerning the ethical use of human volunteers were followed during this research. The study was approved by the Northern Ireland Research Ethics Committee and all participants gave informed consent. Patients were selected from a previously reported cohort [START_REF] Clark | Development of a Diagnostic Genetic Test for Simplex and Autosomal Recessive Retinitis Pigmentosa[END_REF] and all underwent a comprehensive ocular evaluation including Snellen visual acuity testing, visual fields (Humphrey 24-2; Carl Zeiss Meditec, Inc., Dublin, CA, USA), fundal examination, and electro-retinography to establish the diagnosis of RP based on the presence of the typical fundal features (figure 1), visual field constriction and an attenuated or abolished electro-retinogram. DNA was extracted from whole blood using standard protocols (Wizard DNA Purification Kit; Promega, Southampton, UK). 360 DNA samples from the general Northern Ireland population (unscreened for ocular disease) were used as controls.

Sequence capture and next generation sequencing

A custom sequence capture array (Roche NimbleGen, Madison, WI, USA) was designed to target all exons and 100bp of flanking sequence from all genes in which mutations were known to cause RP (Retnet: Retinal Information Network) and selected genes associated with Leber congenital amaurosis, a disease which has phenotypic overlap with RP (table 1). Sequence variants present in greater than 30% of reads at positions covered by at least five reads were classified as novel or known non-coding, synonymous or non-synonymous changes by comparison with the annotated reference sequences. These thresholds were chosen to maximize the detection of real variants whilst minimising artefactual calls and were informed by observation of the detection rate of known SNPs. Identified sequence variants were annotated according to the guidelines published by the Human Genome Variation Society.

Determination of pathogenicity

The potential pathogenicity of sequence variants identified by high-throughput DNA sequencing was assessed as follows. Firstly, known SNPs were excluded as described above.

The frequency of non-synonymous sequence variants in 360 population controls (720 chromosomes) was then determined by either Sequenom mass array genotyping (Sequenom, San Diego, CA, USA) or high resolution melt curve analysis on a LightCycler 480 (Roche Diagnostics, Basel, Switzerland) using Precision HRM MasterMix (Primerdesign, Southamption, UK), according to the manufacturers' protocols. The functional effects of nonsynonymous sequence variants that were absent from or present in less than 1% of the control population were predicted using a range of freely-available computational tools. The conservation of the affected amino acid across species was analyzed using the ClustalW the only affected member in their family (simplex RP) and conventional genotyping and resequencing arrays had failed to detect mutations [START_REF] Clark | Development of a Diagnostic Genetic Test for Simplex and Autosomal Recessive Retinitis Pigmentosa[END_REF] . Such cases present the greatest challenge and these patients were chosen deliberately to assess the ability of targeted sequence capture and high-throughput DNA sequencing to detect novel mutations in the known RP genes.

A custom sequence capture array was designed to capture all exons and 100bp of flanking sequence from a total of 45 genes in which mutations have been reported to cause RP or Leber congenital amaurosis (table 1). When the captured DNA was subjected to highthroughput DNA sequencing approximately 35% of the reads aligned uniquely to the reference sequences, indicative of a greater than 3000 fold enrichment of the targeted sequences. One sample was sequenced in a single flow cell and resulted in high fold coverage of almost the entire targeted region, with a mean coverage of 486 reads and > 20 fold coverage over more than 99% of the target region (figure 2). The four remaining samples were labeled with 4bp tags and pooled prior to sequencing. Although the coverage was reduced to an average of 98 reads, approximately 94% and 88% of the target region was covered by at least 10 or 20 reads respectively (figure 2), demonstrating that sample pooling is a feasible approach to reduce sequencing costs.

Applying a requirement for at least five reads and a minimum variant frequency of 30%, 582 single nucleotide variants were detected, of which 150 were novel variants and 432 were previously reported SNPs. Of the novel variants, 113 occurred in non-coding and 37 in coding sequences (13 of which were non-synonymous) and none affected canonical splice sites. The non-synonymous sequence variants detected included the known homozygous and compound heterozygous mutations, and were all confirmed by Sanger sequencing (table 2).

The known PDE6B homozygous mutation (c.1685G>A; p.G562D) in patient RP148 was detected in all 134 reads which spanned this position (figure 3). For the known CRB1 heterozygous variants 76 of 135 reads (c.2129A>T; p.E710V) and 78 of 130 reads (c.2548G>A; p.G850S) carried the mutated sequence in patient RP179 (figure 3). [START_REF] Bascom | Polymorphisms and rare sequence variants at the ROM1 locus[END_REF] and no mutation was found in the peripherin-RDS gene (MIM *179605).

One novel homozygous missense variant was detected in a simplex patient (RP167; figure 1), in the cyclic nucleotide gated channel beta 1 gene (CNGB1) and was absent in controls (table 2, figure 4). The c.2957A>T sequence variant in CNGB1 exon 29 results in the non-conservative substitution (p.N986I) of the amidic amino acid asparagine (AAC: Asn; N) with the aliphatic isoleucine (ATC: Ile; I). Computational analyses (PolyPhen, PMut and SIFT) predicted pathogenicity (table 2) and evolutionary conservation analysis (ClustalW) showed the N986 position to be invariant across species with CNGB1 orthologs (figure 4). When a recessive RP mutation was detected by high-throughput DNA sequencing the affected patient was counseled about the diagnosis and risk of disease transmission. During this interview the patient revealed that his parents were second cousins, establishing consanguinity in the pedigree. Available family members were screened using Sanger sequencing to assess segregation and the heterozygous state of these family members indicates segregation of the mutation with disease (figure 4).

DISCUSSION

The potential of targeted sequence capture and high-throughput DNA sequencing to improve genetic diagnosis, counseling and deliver personalized medicine is widely acknowledged.

However, applications of this approach as a molecular diagnostic tool have been limited and most commentators report that a number of technical and economic issues must be solved before this technique can be incorporated into routine clinical care. [START_REF] Tucker | Massively parallel sequencing: the next big thing in genetic medicine[END_REF][START_REF] Voelkerding | Next-generation sequencing: from basic research to diagnostics[END_REF][START_REF] Chou | DNA sequence capture and enrichment by microarray followed by next-generation sequencing for targeted resequencing: neurofibromatosis type 1 gene as a model[END_REF][START_REF] Bosch | Keeping up with the next generation: massively parallel sequencing in clinical diagnostics[END_REF] The detection of known and novel RP mutations in this study now establishes DNA pooling and high-throughput DNA sequencing as an effective diagnostic tool for heterogeneous genetic diseases.

Over 40 genes have been identified in RP and there are insufficient phenotypic differences between patients (figure 1) to indicate which gene is likely mutated and guide molecular diagnostics. Inheritance pattern can direct genetic testing, particularly for X-linked RP in which mutations in RPGR (MIM *312610) and RP2 (MIM *300757) predominate, but the majority of patients have no family history. [START_REF] Fishman | Retinitis pigmentosa. Genetic percentages[END_REF] Current molecular testing strategies focus on small sets of genes and apply conventional laboratory techniques based on the relative frequency of mutations and cost-effectiveness of testing. Commercially available genotyping arrays (Asper Ophthalmics, Tartu, Estonia) will only screen for known mutations based on inheritance pattern. [START_REF] Koenekoop | Genetic testing for retinal dystrophies and dysfunctions: benefits, dilemmas and solutions[END_REF] ) Our results demonstrate that high-throughput DNA sequencing can detect previously reported mutations in known RP genes regardless of family history or inheritance pattern.

In three simplex RP patients, conventional genotyping and resequencing arrays [START_REF] Koenekoop | Genetic testing for retinal dystrophies and dysfunctions: benefits, dilemmas and solutions[END_REF][START_REF] Clark | Development of a Diagnostic Genetic Test for Simplex and Autosomal Recessive Retinitis Pigmentosa[END_REF] had failed to detect mutations and high-throughput DNA sequencing was used to look for novel mutations. A novel homozygous mutation was detected in CNGB1 (c.2957A>T) in a patient presumed initially to have simplex RP (figure 1), but who subsequently reported consanguinity in the family (figure 4). The likely functional effects of this novel homozygous variant (figure 4), absence in controls, segregation with disease and consanguinity in the family support its designation as the causative mutation in this patient. CNGB1, which encodes the β-subunit of the rod cGMP-gated channel, is a rare cause of recessive RP (RP45; MIM*600724) with only two mutations previously reported. [START_REF] Bareil | Segregation of a mutation in CNGB1 encoding the beta-subunit of the rod cGMP-gated channel in a family with autosomal recessive retinitis pigmentosa[END_REF][START_REF] Kondo | A homozygosity based search for mutations in patients with autosomal recessive retinitis pigmentosa, using microsatellite markers[END_REF][START_REF] Becirovic | The retinitis pigmentosa mutation c.3444+1G>A in CNGB1 results in skipping of exon 32[END_REF] This mutation would not have been identified by conventional testing strategies because CNGB1 would have been excluded due to the extremely low frequency of mutations.

Therefore, high-throughput DNA sequencing can detect novel variants in addition to any known mutations if applied in previously uncharacterized RP cases. In addition to the immediate benefits for genetic counseling of RP patients high-throughput DNA sequencing provides the opportunity to identify candidates for gene therapy trials. [START_REF] Cremers | Promises and challenges of genetic therapy for blindness[END_REF][START_REF] Smith | Prospects for retinal gene replacement therapy[END_REF] A diagnostic molecular test must deliver a reliable, cost-effective mutational screen which can be performed in a clinical laboratory. Perhaps the most important factor to ensure detection of variants is sufficient coverage. A similar approach for screening of genes involved in RP is being developed by Daiger SP et al [START_REF] Daiger | Targeted high-throughput DNA sequencing for gene discovery in retinitis pigmentosa[END_REF] to target complete genes and additional candidate genes. By focusing on the exonic regions of known RP genes our target region was kept to 359kb, which although well below the capacity of the capture array, meant that coverage was maximized. The minimal coverage of 15-fold recommended for diagnostic implementation of this technology by Hoischen el al [START_REF] Hoischen | Massively parallel sequencing of ataxia genes after array-based enrichment[END_REF] was achieved in 99% of the target region. We demonstrated that sufficient coverage (>15 fold in 90% of the target region) was maintained when four samples, each labeled with a 4 base pair tag or 'bar code' were pooled and sequenced together. Additional cost savings could be achieved by multiplexing of more samples in a first pass screen and if any did not meet specified coverage criteria or no pathogenic variant was detected, further sequencing could be performed to increase coverage. Whilst high-throughput DNA sequencing might have to be outsourced by smaller diagnostic laboratories the data analysis process performed by our group could be readily adopted by any laboratory.

The sheer length of sequence interrogated by targeted capture and massively parallel sequencing inevitably leads to the detection of many sequence variants and the assessment of which of these are potentially pathogenic is a significant challenge. In this study we have focused on the 'low hanging fruit' represented by missense (or nonsense) variants and assessed pathogenicity from predicted functional effects, frequency in controls and segregation with disease. The successful detection of known and novel mutations justifies the introduction of a mutational screen using this approach. More sophisticated analyses are possible, although if a novel non-coding or synonymous variant does not have an obvious effect upon splicing then it is difficult to assess its likely functional effect. It would not be feasible to test all new variants in a control population as part of a routine diagnostic service, however the availability of an increasing number of genome sequences (1000 Genomes Project) will make it easier to classify common and rare variants.

The power of high-throughput DNA sequencing is exemplified by studies in which targeted sequencing of all coding regions, termed the 'exome', successfully identified mutations in two rare Mendelian disorders [START_REF] Ng | Targeted capture and massively parallel sequencing of 12 human exomes[END_REF][START_REF] Ng | Exome sequencing identifies the cause of a mendelian disorder[END_REF] and whole-genome sequencing detected the genetic basis of Charcot-Marie-Tooth disease. [START_REF] Lupski | Whole-genome sequencing in a patient with Charcot-Marie-Tooth neuropathy[END_REF] Unfortunately, the cost of achieving sufficient coverage and the challenges of data analysis mean that the application of exome and whole-genome sequencing as a clinical diagnostic tool is remote. However, we have demonstrated that targeted sequence capture of a disease-specific cohort and subsequent high-throughput DNA sequencing can be used as a cost-effective genetic diagnostic tool.

Figure 1 .

 1 Figure 1. Fundal images and pedigrees of the five retinitis pigmentosa patients. The fundal

Figure 2 .

 2 Figure 2. Sequence coverage of targeted regions. The graph shows the 'completeness' of total

Figure 3 .

 3 Figure 3. Detection of previously identified mutations in two retinitis pigmentosa patients.

Figure 4 .

 4 Figure 4. Detection of a novel mutation (c.2957A>T; p.N986I) in cyclic nucleotide gated

Table 1 . Genes implicated in RP (or Leber congenital amaurosis which shows phenotypic overlap) arranged by inheritance pattern. a

 1 Genes discovered after the design of the sequence capture array (November 2008) and therefore not included.

	Autosomal dominant	Autosomal recessive	X-linked
	CA4	RP9	ABCA4	NRL	RPGR
	CRX	SEMA4A	AIPL1	PDE6A	RP2
	FSCN2	SNRNP200 a	C2ORF71 a	PDE6B
	GUCA1B	TOPORS	CDHR1	PRCD
	IMPDH1		CEP290	PROM1
	KLHL7 a		CERKL	RBP3 a
	NR2E3		CNGA1	RD3
	PRPF3		CNGB1	RGR
	PRPF8		CRB1	RLBP1
	PRPF31		EYS a	RPE65
	PRPH2		GUCY2D	RPGRIP1
	RDH12		IDH3B	SAG
	RHO		LCA5	SPATA7 a
	ROM1		LRAT	TULP1
	RP1		MERTK	USH2A

Table 2 . Classification of non-synonymous sequence variants detected by massively parallel DNA sequencing.

 2 Thirteen non-synonymous variants were detected by massively

	parallel DNA sequencing including 4 mutations (bold text): one known homozygous
	mutation in PDE6B (patient RP148; p.G562D), 2 compound heterozygous mutations in
	CRB1 (patient RP179; p.E710V and p.G850S) and a novel homozygous mutation in CNGB1
	(RP167; p.N986I). a The frequency of non-synonymous sequence variants was assessed in
	360 population controls (720 chromosomes). b Two variants (ROM1: p.R229H; RP1:
	p.L1417P) were classified as rare polymorphisms as they were seen in <1% of controls but
	did not segregate.

c Variants present in at least 1% of the control population were classified as SNPs. d PolyPhen classification: Prob D (probably damaging: high confidence change affects protein function or structure); Poss D (possibly damaging: may affect protein structure or function); and B (benign: lacks phenotypic effect). e PMut classification: N (neutral) and P (pathogenic). f SIFT classification: T (tolerated) and D (damaging).

c.3123C>A; p.H1041Q) and because no second potentially pathogenic variant could be found in the respective genes, they were considered unclassified variants. The non-synonymous variant in ROM1 (MIM *180721; c.686G>A; p.R229H) was found in <1% in controls but was a previously reported polymorphism

ACKNOWLEDGEMENTS

The authors thank Clive Wolsley for his assistance with retinal electrophysiology, Gareth McKay for assistance with control DNA samples and Justin O'Neill for his work in recruiting the families. This research was supported by the Health and Social Care Northern Ireland R&D of the Public Health Agency (Project grant RRG 4•43).

FUNDING

This research was supported by the Health and Social Care Northern Ireland R&D of the Public Health Agency (Project grant RRG 4•43).

COMPETING INTERESTS

The authors have no competing interests

COPYRIGHT LICENCE STATEMENT

The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence (or non-exclusive for government employees) on a worldwide basis to the BMJ Group and co-owners or contracting owning societies (where published by the BMJ Group on their behalf), and its Licensees to permit this article (if accepted) to be published in The Journal of Medical Genetics and any other BMJ Group products and to exploit all subsidiary rights, as set out in our licence. according to hydrophobicity, reveals the high level of conservation of this region. The asparagine (N986) is an invariant residue in a key structural domain: the cyclic nucleotide-binding domain (CNBD). CNBDs are present in cAMP-and cGMP-dependent protein kinases (cAPK and cGPK) and vertebrate cyclic nucleotide-gated ion-channels. Closure of cyclic nucleotide-gated channels occurs when the concentration of cGMP decreases in response to the G-protein-coupled cascade initiated by light. The decrease in cGMP is thus transduced into membrane hyperpolarization.

Previously a homozygous missense mutation in exon 30 of the CNGB1 gene (c. 2978G>T; p.G993V) was identified as a rare cause of recessive RP in a consanguineous French family [START_REF] Kondo | A homozygosity based search for mutations in patients with autosomal recessive retinitis pigmentosa, using microsatellite markers[END_REF] : the residue is indicated (V). Both of these missense mutations in CNGB1 result in amino acid substitutions (p.G993V; p.N986I) in the CNBD and may affect the binding of cGMP to CNGB1.

Modulation of the sensitivity of CNGB1-containing channels is potentially important for light or dark adaptation and can be achieved by phosphorylation of a tyrosine residue (Y983, equivalent to Y1097 in bovine CNGB1 [START_REF] Molokanova | Subunit contributions to phosphorylation-dependent modulation of bovine rod cyclic nucleotide-gated channels[END_REF] ). The p.N986I mutation is only three residues from this tyrosine (P) and is likely to disrupt phosphorylation by altering the phosphorylase binding motif or changing the local 3D protein structure. The loss of 170 amino acids of the distal C-terminus caused by the only other reported mutation in CNGB1 (c.3444 + 1G-->A) [START_REF] Kondo | A homozygosity based search for mutations in patients with autosomal recessive retinitis pigmentosa, using microsatellite markers[END_REF][START_REF] Becirovic | The retinitis pigmentosa mutation c.3444+1G>A in CNGB1 results in skipping of exon 32[END_REF] is indicated by a black line.

Licence for Publication

The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence (or non exclusive for government employees) on a worldwide basis to the BMJ Publishing Group Ltd to permit this article (if accepted) to be published in Journal of Medical Genetics and any other BMJPGL products and sublicences such use and exploit all subsidiary rights, as set out in our licence (http://group.bmj.com/products/journals/instructions-for-authors/licenceforms)


Eric De Silva

J Krishnan
email: krishnan@icex.imperial.ac.uk

Russell Betney

Ian Stansfield

A mathematical modelling framework for elucidating the role of feedback control in translation termination

Keywords: Feedback control, complex control system, translation termination, dynamic analysis, control structure

Translation is the final stage of gene expression where messenger RNA is used as a template for protein polymerisation from appropriate amino acids. Release of the completed protein requires a release factor protein acting at the termination/stop codon to liberate it. In this paper we focus on a complex feedback control mechanism involved in the translation and synthesis of release factor proteins, which has been observed in different systems. These release factor proteins are involved in the termination stage of their own translation. Further, mutations in the release factor gene can result in a premature stop codon. In this case translation can result either in early termination and the production of a truncated protein or readthrough of the premature stop codon and production of the complete release factor protein. Thus during translation of the release factor mRNA containing a premature stop codon, the full length protein negatively regulates its production by its action on a premature stop codon, while positively regulating its production by its action on the regular stop codon. This paper develops a mathematical modelling framework to investigate this complex feedback control system involved in translation. A series of models is established to carefully investigate the role of individual mechanisms and how they work together. The steady state and dynamic behaviour of the resulting models are examined both analytically and numerically.

A c c e p t e d m a n u s c r i p t 1 Introduction

Cell survival and fitness is dependent upon the ability of the cell to alter its complement of enzymes and structural proteins in response to environmental and developmental change. Protein synthesis begins with the transcription of a messenger RNA (mRNA) from a given gene. This mRNA is then translated by the cellular ribosomes to assemble the mature protein. Translation is a complex biopolymerisation process that can be divided into three stages; initiation, in which ribosomes join the mRNA and locate the first, AUG initiation codon; elongation, in which the ribosome moves in a codon by codon stepwise manner along the mRNA, using tRNA carrier molecules to incorporate the appropriate encoded amino acid at each codon, and termination, in which the last, stop codon is recognized by a protein release factor to direct release of the completed polypeptide from the ribosome. In eukaryote organisms, all three stop codons (UAA, UAG and UGA) are recognized by an essential release factor protein, eRF1 [START_REF] Frolova | A highly conserved eukaryotic protein family possessing properties of polypeptide chain release factor[END_REF], which binds ribosomes paused at a stop codon to release the nascent polypeptide. This factor is responsible for the recognition of all in-frame stop codons. Thus the introduction of an early stop codon into an mRNA, even if it lies within the polypeptide encoding portion (through mutation for instance), will cause termination and the production of a truncated protein. In humans such genetic lesions and production of truncated proteins can cause diseases such as cystic fibrosis and Duchenne Muscular Dystrophy [START_REF] Linde | I ntroducing sense into nonsense in treatments of human genetic diseases[END_REF]. However, despite their deleterious nature, premature stop codons within the essential yeast gene encoding eRF1 are not lethal (Stansfield et al 1996, Moskalenko et al., 2003). In this case, an autogenous negative feedback loop controlling eRF1 production is established by the premature stop codon. The early stop codon within the eRF1 gene limits the amount of release factor synthesis by the ribosome. The consequent reduction in eRF1 abundance in turn allows certain tRNA species to compete with the release factor for stop codon recognition. Such tRNAs can either be mutants altered to decode the stop codon (Stansfield et al., 1996) or wild-type, in which case a normal tRNA misreads the stop codon in a translational error process (Moskalenko et al., 2003). Stop codon 'readthrough' or bypass results from tRNA recognition of the stop codon, causing synthesis of full-length release factor (Figure 1). The negative feedback loop thus established plays an important role in maintaining a balance of truncated and full-length release factor synthesis, and is dependent upon the fact that the release factor is recognising a premature stop codon on its own mRNA.

A very similar feedback structure is seen in translation of the mRNA encoding the release factor RF2 in the bacterium Escherechia coli [START_REF] Craigen | Bacterial peptide chain release factos: conserved primary structure and possible frameshift regulation of release factor 2[END_REF][START_REF] Craigen | Expression of peptide chain release factor 2 requires high-efficiency frameshift[END_REF]. The mRNA for release RF2 protein levels in the cell are sufficient to catalyse efficient termination, translation will terminate at this premature UGA stop codon, with a truncated non-functional protein produced. When release factor levels drop, readthrough of the premature stop-codon and frameshifting is triggered, with the consequent production of more full-length RF2 release factor. Here too translation of the mRNA may lead to termination at the premature stop codon, leading to a truncated protein (which will not be able to function effectively as a release factor); alternatively readthrough may occur at the premature stop codon (associated with a frameshift) and further translation will lead to the synthesis of the full length protein. The release factor protein is involved in its own translation termination and also in termination at the premature stop codon. Here again, we have a similar situation with an interplay between feedback at translation termination and the presence of a premature stop codon. Our modelling work in this paper is aimed at elucidating the essential effects of such feedback, the presence of a premature stop codon and their interplay in translation.

For an understanding of the quantitative behavior of such a system it is necessary to build a model to understand the dynamics of such feedback control in translation. Kinetic models of biopolymerization on nucleic acid templates were developed early on [START_REF] Macdonald | Kinetics of biopolymerization on nucleic acid templates[END_REF], [START_REF] Macdonald | Concerning the kinetics of polypeptide synthesis on polyribosomes[END_REF].

These lead to more sophisticated treatments combining experimental results [START_REF] Bergmann | A kinetic model of protein synthesis. Application to hemoglobin synthesis and translational control[END_REF], [START_REF] Heinrich | Mathematical modelling of translation of mRNA in eukaryotes; steady state, time-dependent processes and application to reticulocytes[END_REF]. Further models which included initiation, a lumped elongation and termination [START_REF] Mehra | Insights into the relation between mRNA and protein expression patterns: I. Theoretical considerations[END_REF] have been formulated and extended to codon-based descriptions [START_REF] Mehra | An algorithmic framework for genome-wide modeling and analysis of translation networks[END_REF]. Most recently, a mechanistic model of translation in E. coli has been built to investigate protein synthesis rates (Zouridis and Hatzimanikatis, 2007), (Zouridis and Hatzimanikatis, 2008). Other recent modelling efforts include (Skjandal-Bar et al, 2007) and a detailed study focussing on translation initiation in yeast [START_REF] Dimelow | Control of translation initiation: a model-based analysis from limited experimental data[END_REF]. While these studies provide various insights into the translation process, no modelling work addresses feedback effects in translation termination of the kind discussed above. In this work, the translation termination mechanism and associated feedback effects are the main focus. We develop a modelling framework to study the termination mechanism involving the feedback mechanisms outlined above, and its competition with tRNA-mediated stop codon readthrough. By systematically developing a series of models of increasing complexity, the release factor protein concentration dynamics as well as the role of the feedback In the next section we outline the components and structure of the models in our modelling framework, and formulate the mathematical description. In the following section, we analyze the resulting models in detail both analytically and numerically. In particular we examine carefully the dynamic changes introduced from both the feedback control mechanism, the presence of a premature stop codon and their interplay. Taken together this provides the first steps for a systematic elucidation of such feedback mechanisms and their interplay with premature stop codons in translation.

Method

Our aim is to analyse the release factor dependent feedback regulation by building a series of models that act as a foundation for a systematic analysis of release factor based feedback loop behaviour and its interaction with a premature stop codon, represented in the final model. In model 1 we concentrate on the basic translation and synthesis of the release factor protein. In model 2 we introduce feedback at termination such that complete release factor proteins are required for the termination of the translation of incomplete release factor proteins. In model 3 we examine a variant of model 1 where we introduce a premature stop codon on the release factor mRNA. Finally, in model 4 we incorporate feedback effects mediated by the release factor protein along with the premature stop codon. In this manner we are able to carefully explore the properties of each additional element in turn before combining them. This paper details the model-building process and analyses system dynamics, and steady state behaviour analytically and numerically.

While our work was initially motivated by the specific example of yeast, we formulate and develop a modelling framework which is more general and is also applicable for other cases involving translation termination feedback and premature stop codons in translation. We begin by expanding our description of translation in eukaryotes (see Figure 1). During translation initiation, ribosomes are constantly fed onto the

A c c e p t e d m a n u s c r i p t

mRNA at the 5 end, creating polysomes (mRNAs carrying multiple actively-translating ribosomes). The initiation process begins with the assembly of a so-called Multi-Factor Complex, an assembly of small ribosomal subunits and multiple translation initiation factors. Initiation is complete once the first AUG codon has been recognized on the mRNA (Figure 1). This is followed by a series of elongation steps where the ribosome progresses along the mRNA, one codon at a time, adding a new amino acid at each step to the growing peptide chain. Elongation requires additional elongation factors, as well as a population of transfer RNAs to deliver amino acids in a codon-directed manner. The elongation process is repeated until the ribosome encounters a stop codon (either UAA,UGA or UAG). The stop codon is recognised by a eukaryotic release factor (eRF1) protein acting in concert with a second release factor component, eRF3, allowing the completed protein to detach from the ribosome (Stansfield et al., 1995;Zhouravleva et al., 1995). The complexity of the overall translation system is such that in practice there are many species of mRNA being translated at any given time, competing for the same translational resources (ribosomes, translation factors and tRNAs, amino acids etc). These aspects of the translation are similar in prokaryotes like E.coli, where the release factor RF2 recognizes stop codons to terminate translation.

Our modelling approach is dictated by the need to understand the key features of translation termination, feedback and the effect of introducing a premature stop codon and their interplay. The translation process is therefore modelled in a realistic but simplified way, incorporating detailed descriptions of the biochemistry only when necessary to describe the salient features of the feedback regulation being considered in this work. To understand the essential features of the above factors, it is not necessary to include a detailed explicit description of all elongation steps in translation, as well as detailed descriptions of the biochemistry involved in the initiation step. Therefore these steps are described in a simplified way. Thus rather than describe the elongation phase of translation in a codon-by-codon manner [START_REF] Mehra | An algorithmic framework for genome-wide modeling and analysis of translation networks[END_REF], which describes the translation process in a more detailed manner, we described the elongation part in a lumped fashion. This lumped description will capture the main essential features (steady state and dynamic) of such a detailed model, because the mathematical structure of both descriptions are very similar. It is worth mentioning that our lumped model can be obtained systematically as a reduction of a detailed codon-based model. Any extra biochemical detail in the biochemistry in the initiation and termination steps will not be explicitly described here.

In order to investigate the control structure of the feedback regulation that results from the introduction of a premature stop codon into the release factor mRNA, a deterministic modelling framework has been developed which captures these essential features and this allows for a transparent exploration of A c c e p t e d m a n u s c r i p t the system's attributes. This will act as a platform for understanding the complex control structure and the subsequent development of a stochastic codon-based description involving various mechanistic details.

Through modeling we are able to mathematically describe the autogenous regulation of release factor synthesis, and examine a complex representation of a translation control system in such a way that the emergent system properties are revealed. To this end we describe four primary models of increasing complexity that describe release factor protein translation in the next section.

Model 1

A basic model of translation is illustrated in figure 2 indicating the model variables and parameters involved. This model describes the ribosome joining a free start codon site to result in the occupation of the initiation site, the movement of ribosomes into the polysome and finally to the occupation of the termination site. This final reaction then determines the amount of release factor protein (P ) product produced. Every reaction is modelled using mass-action kinetics.

The mathematical model resulting from the description of the kinetic network above is (see Fig. 2 for details)

dM o dt = (k A .r.E) -(k I .M o ) (1) d(RP ) dt = (k I .M o ) -(k E .(RP )) ( 2
)
dM t dt = (k E .(RP )) -(k T .M t ) (3) dP dt = (k T .M t ) -(k d .P ) (4) dr dt = -(k A .r.E) + (k T .M t ) (5)
In the above RP corresponds to concentration of ribosomes in the polysome pool and M t to the concentration of ribosomes at the termination site. Clearly only the ribosomes at the termination site (which also equals concentration of mRNAs with occupied termination sites) play a role in producing proteins.

Two conservation laws directly result from the model equations. The first captures the fact that the total number of mRNAs (for the release factor protein) is constant and equal to the sum of the concentration of the occupied initial/start codon site (M o ) and the concentration of the unoccupied initial/start codon site E (the equation for E is not written explicitly above). Note that M 0 is the concentration of ribosomes at the initial/start codon site. The second that the total number of ribosomes in the system is conserved. This means that at any time the sum of the concentrations of free ribosomes and those at the start and termination sites and in the polysome pool(s) is a constant, equal to the initial value. This reflects that fact that this model does not describe the creation and breakdown of the ribosomal subunits themselves.

Our goal in this paper is to explore the main qualitative features (both steady-state and dynamic)

arising from the feedback regulation and the presence of a premature stop codon and their interplay, and not on data fitting. We choose a reasonable set of parameters based on estimates and available information.

The value for k A is estimated from the literature (Skjandal- Bar et al , 2007). The rate constant k I is approximately that obtained from a ribosomal scanning rate of ten nucleotides per second [START_REF] Berthelot | Dynamics and processivity of 40S ribosome scanning on mRNA in yeast[END_REF]. The elongation rate, k E , is assumed to be constant and context independent (not influenced by particular sequences of codons). Note that k E implicitly includes the number of codons involved in the elongation phase (k E is essentially inversely proportional to the number of codons, n). The number of codons involved is taken to be 438. An estimate of the degradation rate k d is based on the half-life of the release factor protein measured [START_REF] Belle | Quantification of protein half-lives in the budding yeast proteome[END_REF].

Variables have zero concentration at the beginning of the simulation with the exception of the number of free ribosomes, r 0 and protein, P 0 (note that all mRNAs have unoccupied initiation sites and so

the variable E at time 0 equals the total mRNA concentration). Estimates of total ribosomal numbers per cell vary [START_REF] Beyer | Post-transcriptional expression regulation in the yeast Saccharomyces cerevisiae on a genomic scale[END_REF] but we assume a fraction of this total are available for translation and this number is represented by the number of ribosomes at time zero of the simulations, r 0 . We use estimates of release factor concentrations for estimating initial conditions of the protein P 0 [START_REF] Ghaemmaghami | Global analysis of protein expression in yeast[END_REF]: in this case we estimate the initial concentration of protein to be apprioximately half that number, the reduced number accounting for the fact that these proteins are also involved in liberation of other synthesized protein species. In our simulations, we consider a ribosome to mRNA ratio of 10.

Model 2

In this model feedback is introduced so that the release factor (P ) is needed for the termination and hence synthesis and production of new release factor (P ). We do not include any redundant/extra termination factor here (see conclusions). Figure 3 shows the components of this model.

A c c e p t e d m a n u s c r i p t

The model equations are

dM o dt = (k A .r.E) -(k I .M o ) (6) d(RP ) dt = (k I .M o ) -(k E .(RP )) (7) dM t dt = (k E .(RP )) -(k T .M t .P ) (8) dP dt = (k T .M t .P ) -(k d .P ) (9) dr dt = -(k A .r.E) + (k T .M t .P ) (10)
The equations above reflect the presence of the protein in translation termination results as seen in equations ( 8), ( 9) and ( 10) above. This also means that in the above equations the termination rate constant k T is a second order rate constant. It is also worth noting that with feedback, at the termination site M t , the termination involves an interaction of release-factor protein with the ribosome, resulting in the release of the newly synthesized release factor protein and then being released itself.

Model 3

The introduction of a premature stop codon into the release factor mRNA means that the elongation stage can be divided into two stages, which is reflected in the two polysome pools (RP 1 ) and (RP 2 ). The protein being synthesized may be released at the first termination site (giving rise to the truncated protein P t ), or following readthrough of the premature stop codon, may continue elongation and eventually be released at the regular termination site. The readthrough of the premature stop codon is mediated by a tRNA, denoted by T. In the case of yeast, a mutated serine tRNA was responsible for this readthrough (Stansfield et al 1996). In the case of E.coli, the readthrough of the (naturally occurring) premature stop codon is associated with a frameshift. The schematic for this model is illustrated in figure 4.

The time variation of various concentrations is obtained by simulating the above equations. T is a parameter whose concentration is assumed to constant and k t1 is a second order rate constant associated with readthrough. The effect of tRNA concentrations on the dynamics of the system is investigated by simulating for a range of tRNA concentrations. Thus the effect of varying T is that of varying the rate constant of the quasi-first-order reaction involving readthrough at M t1 . The new coupled differential equations describing Model 3 are given as follows:

A c c e p t e d m a n u s c r i p t

The equations depicting this process are given by

dM o dt = (k A .r.E) -(k I .M o ) (11) d(RP 1 ) dt = (k I .M o ) -(k E1 .(RP 1 )) ( 12
)
dM t1 dt = (k E1 .(RP 1 )) -(k T a .M t1 ) -(k t1 .M t1 .T ) ( 13
)
dP t dt = (k T a .M t1 ) -(k dt .P t ) (14) d(RP 2 ) dt = (k t1 .M t1 .T ) -(k E2 .(RP 2 )) ( 15
)
dM t2 dt = (k E2 .(RP 2 )) -(k T b .M t2 ) (16
)

dP dt = (k T b .M t2 ) -(k d .P ) ( 17
)
dr dt = -(k A .r.E) + (k T a .M t1 ) + (k T b .M t2 ) (18)
These equations describe the process depicted above. Note that protein synthesis can terminate either at the premature stop codon or the regular stop codon, associated with which are the termination rate constants k T a and k T b respectively. Note that again, as before the elongation constants k E1 , k E2 implicitly include the number of codons in each elongation phase (inversely proportional to the number of codons).

For the purposes of simulations, we assumed a premature stop codon at the 100th codon location. In the simulations, for specificity we assumed that the two termination rate constants were equal and that the degradation rate of the truncated protein was less than that of the full length protein.

Model 4

This model includes both a premature stop codon as in Model 3 (and the associated tRNA which mediates readthrough) and feedback effect on translation termination as in Model 2 and is illustrated in figure 5.

Note that the full length protein is responsible for termination both at the regular as well as the premature stop codons.

The feedback results in an extra P term at both termination sites seen in equations ( 21), ( 22), ( 24), (25), and (26) below. As before the product of k t1 and the parameter T represents a pseudo-first order rate constant.

A c c e p t e d m a n u s c r i p t

dM o dt = (k A .r.E) -(k I .M o ) (19) d(RP 1 ) dt = (k I .M o ) -(k E1 .(RP 1 )) ( 20
)
dM t1 dt = (k E1 .(RP 1 )) -(k T a .M t1 .P ) -(k t1 .M t1 .T ) ( 21
)
dP t dt = (k T a .M t1 .P ) -(k dt .P t ) (22) d(RP 2 ) dt = (k t1 .M t1 .T ) -(k E2 .(RP 2 )) ( 23
)
dM t2 dt = (k E2 .(RP 2 )) -(k T b .M t2 .P ) ( 24
)
dP dt = (k T b .M t2 .P ) -(k d .P ) ( 25
)
dr dt = -(k A .r.E) + (k T a .M t1 .P ) + (k T b .M t2 .P ) (26)
Numerical integration of the model equations described above allows exploration of their dynamical behaviour. Simulations are carried out on MATLAB 1 using the ode15s stiff solver. A sample set of parameters are chosen to describe translation. While an exploration of a few important parameters has been carried out, a more exhaustive parametric analysis will be examined at a later time. The dynamical behaviour of the model is also investigated via examination of the relevant analytical expressions. These approaches are supplemented by the use of chemical reaction network theory [START_REF] Feinberg | Chemical reaction network structure and the stability of complex isothermal reactorsI. The deficiency zero and deficiency one theorems[END_REF] through the Network Analyst software part of the Chemical Reaction Network Toolbox package 2 to investigate whether the differential equations that describe the model allow oscillations and/or multiple steady states, or whether such scenarios are ruled out based on the structure of the relevant reaction networks. In certain cases chemical reaction network theory is able to provide additional information about the possible asymptotic dynamic states of the network, based on the structure of the network alone (valid for reaction networks with mass action kinetics).

Results

In this section we present results from the analysis of the various models. Numerical simulations are performed for each of the dynamical models. The effects of some important parameters are also examined.

A sequential examination of models allows for a clear analysis of the features of each model, and hence

Model 1

Some results from the simulations are shown in Figure 6 which shows how the concentrations of the occupied ribosome pool (RP ) and protein (P ) vary with time within the Model 1 framework. The system reaches a steady state and simulations starting with different initial conditions provide no indication of other steady states (see below). As part of the simulations, the effect of some parameters was examined on the steady states. It was found, for instance, that the protein concentration P increased as the termination rate constant k T increased. For example, protein production P showed a weak increase as k T is increased by a factor of 10 (keeping all other parameters fixed). The protein concentration was found to be comparitively more sensitive to a ten-fold increase in total ribosome concentation.

While numerical sensitivity analysis reveals some aspects of a parameter's influence, one parameter at a time, greater insight can be obtained by obtaining analytical expressions for the steady state(s) of the model equations. In order to obtain the steady state concentrations of different variables the right hand side of equations ( 1)-( 5) are set to zero. This results in the following equation for (P ):

P = k T .M t k d ( 27
)
where

M t = k E .(RP ) k T = k A .r.E k T (28)
The above equations along with the conservation conditions for the number of free ribosomes and the total number of mRNAs result in:

M o + (RP ) + M t + r = K (29)
where K is a constant equal to the number of free ribosomes initially, and:

E = (MR) -M o ( 30
)
where MR is the total concentration of mRNAs (a constant). Then from ( 27) and ( 29) we can write P :

P = k T k d K 1 + k T . 1 k I + 1 k E + 1 k A .E (31)

A c c e p t e d m a n u s c r i p t

We now write the variables entirely in terms of P and rate constants:

M t = α 1 .P ; (RP ) = α 2 .P ; M o = α 3 .P ; r.E = α 4 .P (32)
where

α 1 = k d k T ; α 2 = k d k E ; α 3 = k d k I ; α 4 = k d k A (33)
The conservation equations can then be written as:

E = (MR) -α 3 .P and r = K -(α 1 + α 2 + α 3 ) .P (34)
Using r.E = α 4 .P , P can be written as a quadratic equation:

(α 1 + α 2 + α 3 ) .α 3 .P 2 -(K.α 3 + (α 1 + α 2 + α 3 ) .(MR) + α 4 ) .P + K.(MR) = 0 (35)
The solution to the above equation is:

P = β ± β 2 -4. (α 1 + α 2 + α 3 ) .α 3 .K.(MR) 2. (α 1 + α 2 + α 3 ) .α 3 (36)
where

β = (K.α 3 + (α 1 + α 2 + α 3 ) .(MR) + α 4 ) (37)
This expression suggests that there will be two real and positive roots for P in general. Incidentally, for the model 1 parameter values this results in two roots for P , P =0.07 and P =0.01 moles/vol. While the simulations reveal only one steady state it is important to examine the presence of the other steady state. It turns out that the second steady state results in negative values for the concentrations for other model variables, such as those of free initiation sites, E. This is seen directly from the solution for the steady state of E (writing the equation for the steady state of E, by eliminating other variables):

E = -γ ± γ 2 + 4.α 4 .(MR) α 3 . 1 + 1 α3 . (α 1 + α 2 ) 2. 1 + 1 α3 . (α 1 + α 2 ) (38)
where

γ = 1 α 3 (α 4 -α 1 .(MR) -α 2 .(MR)) -K -(MR) (39)
Equation ( 38) immediately reveals that one of the two roots for E (and hence P ) is inadmissible because it is negative. Examination of r likewise reveals a positive and negative root. These results are confirmed by (positive) values the rate constants might have. In summary, the analytical expressions reveal the exact dependence of different steady state concentrations on different parameter values explicitly. This analysis is useful in understanding more complex subsequent models. Furthermore, the analysis directly reveals the presence of a single admissible steady state (all variable concentrations non-negative).

Model 2

Model 2 involves feedback effects in translation termination at the regular termination site, Once again multiple simulations are performed in order to see how the concentrations of different variables vary with time. Figure 7 shows how the concentrations of the occupied ribosome pool (RP ) and protein (P ) vary with time. Simulations were repeated by changing initial conditions. In all cases a single steady state was observed. This is confirmed by the Chemical Reaction Network Toolbox software which under mass-action kinetics calculates a deficiency of order one [START_REF] Feinberg | Chemical reaction network structure and the stability of complex isothermal reactorsI. The deficiency zero and deficiency one theorems[END_REF], concluding that taken with mass action What is apparent from the structure of the Model 2 chemical network is that the protein product P is needed for its own production, so that the initial conditions for P must be non-zero (as is the case for a living cell). Even very low P concentrations at the beginning of the simulations result in the steady state value of P being the same as for much higher initial concentrations. We once again have five equations at steady state which allow us to write P as:

P = k E .(RP ) k T .M t = k A .r.E k T .M t = k I .M o k T .M t ( 40
(K -α 3 .P -α 2 .P -α 1 ) . ((MR) -α 3 .P ) = α 4 .P (42)
producing the quadratic:

(α 2 + α 3 ) .α 3 .P 2 + δ.P + K.(MR) -α 1 .(MR) = 0 (43)
where

δ = α 1 .α 3 -α 4 -(α 2 + α 3 ) .(MR) -α 3 .K (44)
Once again there is the possibility of two real roots:

P = -δ ± δ 2 -4.α 3 .(MR). (α 2 + α 3 ) . (K -α 1 ) 2. (α 2 + α 3 ) .α 3 (45)
So the above expression shows that P has two feasible (positive) roots for our basal parameter values and similar examination of the corresponding solutions to E and r shows that both have one feasible and one infeasible root. Thus one of the roots for P above can be ruled out in general as it will result in infeasible solutions for E and r.

While the feasible root for P represents a stable steady state (for our basal parameter values) there also exists another steady state for when there is no protein (P =0). Here, there is no protein to aid in its own production and from the steady state equations one can see that M o =0 and (RP )=0, and further, no protein is produced. However, M t , is nonzero. That this is the case is also seen when we simulate model 2 with P =0 initially. This is in contrast to the other steady state where P = 0 and M t = k d /k T . By varying parameters in simulations and significantly increasing the protein degradation rate constant k d or significantly lowering the termination constant k T results in a P =0 being observed as the (stable) steady state. So the introduction of feedback in translation termination has the consequence that if the protein degradation rate is too high or the termination rate constant too low, P cannot sustain its own production.

A c c e p t e d m a n u s c r i p t

Thus, in general there are two steady states, one corresponding to zero protein concentration and another non-zero steady state. Under basal conditions the non-zero (positive) steady state is stable while the zero steady state is unstable. Decreasing the termination rate constant sufficiently results in a scenario

where the P = 0 steady state is stable and the non-zero steady state is in fact inadmissible (negative concentrations of entities in the network). The transition between these scenarios occurs via a transcritical bifurcation (discussed in detail both analytically and numerically elsewhere).

In summary, simulations of the model for basal parameter values reveals a non-zero steady state of protein concentration. The presence of an additional steady state corresponding to zero protein concentration is inferred analytically and confirmed by simulations. For the basal parameter values chosen here the non-zero steady state is the only stable steady state.

Model 3

In Model 3 we introduce a premature stop codon (chosen to be the 100th codon site in our simulations).

The model contains the concentration of a tRNA protein (T ) which is responsible for readthrough at the premature stop codon. In bacterial systems, there are similar factors which facillitate readthrough, with an associated frameshift. There is no feedback in Model 3 so the effects of multiple stop codons can be explored separately, independent of feedback. Multiple simulations are run to examine the changes in variable concentrations with time. Figure 8 shows how the concentrations of the occupied ribosome pool before (RP 1 ) and after (RP 2 ) the premature stop codon and the full-length (P ) and truncated protein (P t ) vary with time within the Model 3 framework. The steady state is approached exponentially and there is no indication of other steady states. Chemical reaction network theory is unable make predictions about Model 3 due to the structure of the network (deficiency =2).

We now examine the steady state of this more complex network in detail. In this model the variables of interest are the concentration of full-length, P , and truncated, P t , protein. In Model 3 a particular aspect of interest is the ratio of of full-length to truncated protein which can be obtained as:

P = k T b .M t2 k d and P t = k T a .M t1 k dt (46)
From this and the steady state equations for the model, we find that

P P t = k dt k d . M t2 M t1 = k dt .k t1 .T k d .k T a (47)

A c c e p t e d m a n u s c r i p t

We notice straightaway even without explicitly solving the steady state equations that the ratio of full-length to truncated proteins depends only on rate constants and the tRNA concentration, T .

There are now eight steady state equations derived from eqns. ( 11)-( 18). concentration, T . From these eight steady state equations all the variables can be written in terms of only P and the rate constants: where:

M t2 =
α 1 = k d k T b ; α 2 = k d k E2 ; α 3 = k d k t1 .T ; α 4 = k T a .k d k E1 .k t1 .T + k d k E1 ; ( 50
)
α 5 = k T a .k d k dt .k t1 .T ; α 6 = k T a .k d k I .k t1 .T + k d k I ; α 7 = k T a .k d k A .k t1 .T + k d k A ; (51)
The conservation equation E + M o = (MR) still holds while the conservation rule for the total number of ribosomes is now:

M o + (RP 1 ) + M t1 + (RP 2 ) + M t2 + r = K (52)
Combining these with the expression for r.E (eqn. 48) we can write:

(Kθ.P ) . ((MR)α 6 .P ) = α 7 .P where θ = α 1 + α 2 + α 3 + α 4 + α 6 (53)

From the above we obtain a quadratic expression:

θ.α 6 .P 2 -(K.α 6 + α 7 + θ.(MR)) .P + K.(MR) = 0 (54)
for which the two roots are given by:

P = -κ ± κ 2 -4.θ.α 6 .K.(MR) 2.θ.α 6 where κ = -(K.α 6 + α 7 + θ.(MR)) (55)
For the basal parameters this results in two positive feasible roots for P . Solving in a similar manner for E results in a positive and a negative root. The feasible root is easily identified and agrees with results observed in the simulations. Thus in summary there exists only one admissible steady state, for reasons very similar to model 1.

A c c e p t e d m a n u s c r i p t

Readthrough efficiency

Since a key factor which determines the concentrations of full-length and truncated proteins is the tRNA concentration T (which controls the degree of premature stop codon readthrough) we examine the effect of this parameter in detail. We already know from above that the ratio of P to P t is proportional to T .

However it is important to determine how these variables themselves change with T . To carefully explore the effect of T , including various limiting cases, we have run multiple simulations for varying amounts of T over a wide range and the results are shown in figure 9(a). We find that the steady state concentration of P increases with increasing T eventually saturating at high T . The steady state concentration of P t does the opposite, saturating at low T and decreasing at high T .

In order to examine what the analytical expressions indicate about P and P t dependency upon tRNA concentration (T ) we use equation ( 55) for P and solve for a range of values of T matching those in the simulations. The feasible root for P is easily identified. The steady state for P t is then easily obtained and both P and P t versus T are shown in figure 9(b). The analytical plots of P and P t versus T match those obtained in the simulations. They show the gradual increase in P and then a leveling off for higher T values. Conversely, the P t concentration plateaus at low T and falls off as T increases (falling at the same time as P levels off). These trends can be understood as follows. The P concentration saturates at high T because the flux of ribosomes through the pathway leading to the premature termination site (and beyond the premature termination site) limit the rate of production of protein P at high T . At high T , readthrough dominates the possibility of termination at the premature stop codon. Thus, very little truncated protein is produced, The situation is different at low T . Here the termination at the premature stop codon dominates the possibility of readthrough. Thus, the concentration of truncated protein P t saturates for low T (being limited by other factors), while the concentration of P approaches zero

From the simulated and analytical results we see that the concentrations of full-length and truncated protein equal one another when:

T = k T a .k d k dt .k t1 (56)
In general the asymptotic concentrations for P t (at low T ) is different from the asymptotic concentration for P (for high T ). Note that the degradation rates of the proteins are in general different.

To summarize: when T is very large the rate of protein production is then limited by the speed at which ribosomes can move through the system. Likewise, when T is very small the rate of production of truncated protein is limited by the speed at which ribosomes pass through and terminate at the premature

Model 4

In this model we introduce feedback at termination to the premature stop codon system of model 3. We perform some preliminary analysis to understand the interaction of this feedback mechanism with a premature stop codon. Multiple simulations are run to examine how variable concentrations change with time. Figure 10 shows how the concentrations of the occupied ribosome pool before (RP 1 ) and after (RP 2 ) the premature stop codon and the full-length (P ) and truncated protein (P t ) concentrations vary with time within the Model 4 framework. Steady state is approached exponentially and there is no indication of oscillatory behaviour. Once again chemical reaction network theory is unable to analyse Model 4 due to the structure of the network. Calculating the deficiency of this reaction network (a mass-action kinetic network) reveals that the deficiency is two [START_REF] Feinberg | Chemical reaction network structure and the stability of complex isothermal reactorsI. The deficiency zero and deficiency one theorems[END_REF] and therefore reaction network theory is unable to make any strong predictions.

We have explored the effect of varying the initial concentration of P in the system to see the effect this has on P and P t production under a range of T concentrations. Increasing the initial P by two orders of magnitude and decreasing initial P by two orders of magnitude result in the same P and P t concentrations at steady state.

The large parameter space in this model indicates that examination of the analytical solutions would be beneficial. As in model 3 there are eight steady state equations and we can write the variables in terms of P and other rate constants (assuming P is non-zero) to give:

M t2 = α 1 ; RP 2 = α 2 .P ; M t1 = α 3 .P ; RP 1 = α 4 .P 2 + α 5 .P ( 57
)
P t = α 6 .P 2 ; M o = α 7 .P 2 + α 8 .P ; r.E = α 9 .P 2 + α 10 .P (58)

where

α 1 = k d k T b ; α 2 = k d k E2 ; α 3 = k d k t1 .T ; α 4 = k T a .k d k E1 .k t1 .T (59) α 5 = k d k E1 ; α 6 = k T a .k d k dt .k t1 .T ; α 7 = k T a .k d k I .k t1 .T (60) α 8 = k d k I ; α 9 = k T a .k d k A .k t1 .T ; α 10 = k d k A (61)

A c c e p t e d m a n u s c r i p t

Using the the conservation conditions described previously one can write:

r.E = α 9 .P 2 + α 10 .P (62) and this along with the above equations results in a complicated quartic equation for P , some of whose coefficients depend on T :

Thus, in contrast to model 2 where we derived a quadratic equation for P in model 4 we get a very complex quartic equation. The complexity of this equation means that it cannot be easily solved analytically. However, just as in model 2, in addition to the steady state(s) described above, there is a second steady state which corresponds to zero protein concentraion (P = 0) as is immediately seen from the steady state equation for P , equation ( 25):

k T b .M t2 .P = k d .P (63)
This steady state is also associated with zero concentration of truncated protein P t , since the full length protein is needed to release truncated protein at the premature stop codon. This steady state is unstable for our basal parameter values (as well as other T concentrations examined below). Just as in model 3 we examine the role of T concentration in determining the concentrations of full-length and truncated proteins, this time in the presence of feedback.

Readthrough efficiency

Now we look at the interplay of feedback and readthrough by again simulating the steady state concentrations for a range of tRNA (T ) concentrations. P and P t can be written in terms of one variable, in this case M t1 :

P = k E2 .(RP 2 ) k T b .M t2 = k t1 k d .M t1 .T (64)
and:

P t = k T a k dt .k d . (M t1 ) 2 .k t1 .T (65)
The P/P t ratio can also be written in terms of T using the steady state equations:

P P t = k dt .k t1 k T a .k E2 . 1 (RP 2 )
.T (66)

In figure 11(a) we see that increasing T concentration in our simulations results in increasing P concentration up to a point after which the value saturates again (at steady state) in similar fashion to model

A c c e p t e d m a n u s c r i p t

3 (figure 9). The P t concentration however increases as T increases, peaks, then begins to decrease in concentration at high T levels. This biphasic shape results from expressions ( 59) and ( 61) written out fully below:

P t = k T a .k d .P 2 k dt .k t1 .T ( 67
)
where P t is inversely proportional to T but proportional to P 2 . This is in striking contrast to the corresponding expression for model 3, where P t is proportional to P (and inversely proportional to T ).

This P 2 /T factor and the inherent saturation of P results in the shape we (discussed below). Thus, as

T becomes large P saturates and P t approaches zero, as seen in figure 11(a). This can be understood physically as follows: when T is very high the flux through the readthrough path is dominant with very little premature termination. However, in this case the P concentration saturates as it is limited by various other rate constants in the pathway.

The other asymptotic limit of low T concentration is also worth examining carefully. We find that for small T the P concentration is essentially linearly proportional to T , and hence from equation ( 69), so is P t . Thus for low T concentrations both P and P t concentrations approach zero. This is easily understood physically as follows: when T concentrations are very low the readthrough at the premature stop codon (modulated by T ) is the rate-limiting step in the production of P . However, since P is itself responsible for termination at the premature stop codon also, this results in low truncated protein concentrations as well. Thus the biphasic dependence of P t concentration on tRNA concentration can be traced to the effect of feedback in translation termination. Figure 11(b) shows the way in which full-length protein varies with truncated protein for a range of T values.

Conclusions and Discussion

The goal of the paper was to understand the effect of the control structure involving feedback at translation termination and the interplay with a premature stop codon. Such feedback regulation has been observed in translation in multiple systems. In yeast there is feedback regulation of a release factor on its own translation termination process, and experiments have been conducted to introduce premature stop codons and analyze the effect of natural feedback in regulating translation in the presence of premature stop codons (Stansfield et al., 1996;Moskalenko et al., 2003). Similarly in E.coli, the release factor RF2 contains a naturally occurring premature stop codon and a similar feedback regulatory structure exists [START_REF] Craigen | Bacterial peptide chain release factos: conserved primary structure and possible frameshift regulation of release factor 2[END_REF][START_REF] Craigen | Expression of peptide chain release factor 2 requires high-efficiency frameshift[END_REF]). In the E.coli case, the readthrough of the premature stop codon is

A c c e p t e d m a n u s c r i p t

associated with a frameshift. Our approach is aimed at understanding the essential elements underlying the control mechanism discussed above. To that end, we formulate a series of models, which allows for the careful and systematic consideration of feedback,the effect of a premature stop codon and their interaction.

Each model was investigated computationally and analytically. Dynamic simulations of model 1 along with analytical work (a simple basic model of translation) indicated a single steady state and its dependence on parameters. The second model involves the introduction of feedback at translation termination. In addition to the positive steady state, there exists a second steady state where there is zero protein produced. This steady state was unstable for our basal parameter values. A sufficient reduction in the termination rate constant or increase in protein degradation rate constant could cause the non-zero steady state to be inadmissible (concentrations of some entities negative) and the zero steady state to be stable. A transcritical bifurcation lies at the transition between these two scenarios. In our modelling, we did not incorporate any residual/redundant additional protein acting as another release factor (none have been identified biologically in these systems). The introduction of such a factor would lead to the breaking of the transcritical bifurcation. However, the broad nature of variation of steady state protein concentrations as a function of parameters would remain the same.

The addition of a premature stop codon (to model 1) and the tRNA responsible for readthrough (model 3) introduces two possible pathways that translating ribosomes can take. At low tRNA concentrations, as tRNA concentration decreases the concentration of full-length protein decreases and that of the truncated protein saturates. On the other hand, as the tRNA concentration becomes large, the full-length protein concentration approaches a constant value and the truncated protein concentration decreases. In model 4 (incorporating a premature stop codon and feedback in translation termination) there exists one positive steady state and the zero steady state of full length protein concentration. (which is unstable for our basal parameter values). Here (unlike model 3) the dependence of the truncated protein concentration on tRNA concentration results in a biphasic relationship: at low tRNA concentrations, the truncated protein concentration approaches zero. This directly stems from the fact that full-length protein is required for termination at the premature stop codon. It should be mentioned that in model 4, we have implicitly assumed a premature stop codon location, which would result in a (sufficiently) truncated protein which is unable to function as a release factor itself. It may be noted that in yeast, a truncation of eRF1 (a 438 amino acid protein) by more than 32 amino acids leads to a protein which is unable to function as a release factor and results in cell death [START_REF] Eurwilaichitr | The C-terminus of ERF1 defines a functionally important domain for translation termination Saccharomyoces cerevisiae[END_REF] In essence, our analysis suggests a few important points. The effect of feedback regulation involved in translation of release factors suggests that the degradation of the release factor protein must be small enough to allow for sufficient release factor protein concentrations to be realized. If either due to internal or external factors, the release factor proteins are degraded faster than normal, it can lead to a very substantial depletion of this protein, which very severely limits synthesis of the protein itself. This fact taken together with the experimentally known fact that knock-outs of these release factor proteins prove lethal indicate that these release factor protein degradation rates must be kept with certain bounds for the survival of the cell. This same insight is relevant for translation in the presence of premature stop codons for the release factor. In the presence of a premature stop codon, the feedback effect allows for a better readthrough of the premature stop codon, allowing release factor protein to be produced. The importance of the degradation rate of the release factor protein not being too high has been demonstrated in experimental observations of yeast carrying sup45 alleles with premature stop codons, thus involving the kind of feedback discussed [START_REF] Kiktev | The paradox of viable sup45 STOP mutations: a necessary equilibrium between translation readthrough, activity and stability of the protein[END_REF]. In this work the authors establish the stability of the encoded mutant form of eRF1 as a crucial factor in regulating the negative feedback loop and maintaining the viability of these strains.

In general the above analysis also suggests the possibility of a redundant protein involved in translation termination in these systems. The above analysis provides a quantitative description of the feedback regulation involved (in particular) in the translation of release factors with premature stop codons.

It should be mentioned that another biological mechanism, that of nonsense-mediated decay, operates in eukaryotes to degrade mRNA with premature stop-codons (Ruiz-Echevarria et al., 1996;[START_REF] Cao | Computational modelling and experimental analysis of nonsense-mediated decay in yeast[END_REF]. A mechanistic description of this process and its interaction with the feedback regulation described in this paper is not attempted here. However, previous work shows that the NMD pathway effect does influence the nature of translational feedback loops operating on eRF1 alleles with premature stop codons and that inactivation of NMD increases the viability of such alleles [START_REF] Chabelskaya | Inactivation of NMD increases viability of sup45 nonsense mutants in Saccharomyoces cerevisiae[END_REF].

There are various other aspects to be explored further in our modelling framework. A detailed parametric sensitivity/bifurcation analysis with respect to various parameters will be performed. Further modelling will explore some additional mechanistic detail in understanding this feedback mechanism in yeast. One detail which needs to be incorporated is a more comprehensive mechanistic description of termination involving the protein (eRF1) in a complex with another protein (eRF3). Future modelling work will also involve an explicit description of the dynamics of tRNA concentration. While our modelling has focussed on a particular feedback control structure, similar but more complex feedback structures will be examined subsequently. One example is that of a system very similar to E.coli is represented by the release factor-2 (RF2) gene of the bacterium Salmonella typhymurium, which like E.coli, also has a +1 frameshift premature stop codon lies at the regulatory frameshift site. Here, in addition to this frameshift site premature stop codon, an RF2 mutation was isolated that resulted in the introduction of an additional in-frame premature stop codon downstream of the initial frameshift site termination codon. This is an example of an additional complexity built around the regulatory structure studied here.

There are other biophysical effects which could interact with the feedback regulatory structure discussed here. One example is the recycling of detached ribosomes back to the initiation site. A detailed analysis of this biophysical phenomenon and the role of mRNA loop formation in translation regulation is considered by modelling these effects on an isolated mRNA chain [START_REF] Chou | Ribosome recycling, diffusion, and mRNA loop formation in translational regulation[END_REF]. In the context of our feedback regulatory mechanism, such additional effects could play modulatory roles. For instance if the premature stop codon was present early along the mRNA chain, then this could lead to enhanced recycling and generally higher ribosomal throughput, when compared to the case where the premature stop codon was present near the termination site. Thus in this case one might expect higher concentrations of both full length and truncated proteins. However a careful and systematic analysis needs to be performed.

Another step is to expand the modelling description in this paper into a detailed codon-based model (mean-field). The modelling work in this paper provides non-trivial insights into the effect of the control structure regulating translation termination and its interplay with the premature stop codon. Since the mathematical structure of the models presented here and corresponding detailed codon-based models are very similar, it is expected that the essential insights carry through in that case too. Eventually, this description needs to be expanded to a stochastic codon-based model. Such models can be formulated as Totally Asymmetric Simple Exclusion Process (Shaw et al., 2003;[START_REF] Chou | Clustered bottlenecks in mrna translation and protein synthesis[END_REF]Shaw et al., 2004a). This would allow for a more detailed and realistic description of various aspects of translation (for example multiple codon coverage of ribosomes). Further, additional aspects of ribosomal mechanochemical cycle have been also incorporated in models in the TASEP setting [START_REF] Basu | Traffic of interacting ribosomes: effects of single-machine mechanochemistry on protein synthesis[END_REF].

The role of stochasticity can be carefully and systematically explored in a TASEP setting, in conjunction with other mechanistic details described earlier. The modelling work in this paper provides a basis for systematically assessing the effects of these additional complexities and their interaction with the feedback control structure. It also provides a platform for more detailed mechanistic modelling of individual steps in translation in different systems with similar feedback control structure. In general, we see that higher tRNA concentrations lead to higher concentrations of full length protein and lower concentrations of truncated protein.

  in E.coli indicates the presence of a (naturally occurring) premature stop codon. The RF2 mRNA's open reading frame is split into two parts by a premature UGA stop codon. The reading frame downstream of the stop codon lies in the +1 frame relative to the upstream, 0 frame portion. When

  examined carefully. The resulting mathematical modelling framework is a tool for carefully elucidating the role of such feedback control mechanisms in the context of translation, and their interplay with other processes in translation. While this modelling work was initially motivated by the specific example of dynamics of the release factor protein eRF1 in yeast, we formulate a more general modelling framework which is relevant to translation with feedback regulation in translation termination, premature stop codons or their combination. Thus the methods and insights are relevant for translation in other systems with similar feedback control structures and/or premature stop codons.

  Network Toolbox software which concludes that taken with mass action kinetics, the network cannot admit multiple positive steady states or a degenerate steady state no matter what

  kinetics, the network cannot admit multiple positive steady states no matter what (positive) values the rate constants might have. The simulations show that decreasing the concentration of mRNA results in less protein produced and more free ribosomes. Increasing the mRNA concentration to exceed that of the ribosomes results in little change in protein production (compared to basal values, r/(MR)=10) but a decrease in the number of free ribosomes. As part of the exploration of Model 2 the initial variable values at time zero in Model 2 simulations are populated with the steady state values from Model 1, however this results in no change to the steady state values.

  )Interestingly, for the equation describing the rate of change of P (equation 9):k T .M t .P = k d .P(41) that the concentration of ribosomes at the termination site depends only on the ratio of protein degradation to termination rate constants (M t = k d /k T ). In model 2 the steady state variable concentrations can be solved for analytically just as before. In order to determine the solutions to the equations we once again write the variables in terms of just P and the rate constants (α 1 , α 2 , α 3 and α 4 are the same as in model 1). This set of equations are the same as in model 1 with the exception of the M t variable, and so using the conservation conditions E + M o = (MR) and M o + (RP ) + M t + r = K (which still apply) one can write r.E = α 4 .P (eqn. 32) as:

  These limiting cases provide insight into the role of T in determining the steady state P and P t concentrations.

  factor expression via a negative feedback loop. As in the E.coli RF2 gene, a

Fig. 1 .

 1 Fig. 1. Schematic of the translation process. A schematic of the translation process is shown in (a) indicating the gradual building of the polypeptiode chain and release at the final stage. The presence of a premature stop codon is shown in (b), along with the feedback loop involved in the release factor protein acting to terminate translation at the premature stop codon (the action of the release factor at the regular stop codon is not shown).

Fig. 2 .

 2 Fig. 2. Schematic of model 1 This figure depicts a schematic of model 1, with ribosome binding to a free initiation site, then progressing to the polysome pool, and eventually terminating to produce the synthesized protein (and releasing the ribosome, which is not shown explicitly). Basic Model 1 variables: r = concentration of free ribosomes; M o = concentration of ribosomes at start/initiation codon site; (RP ) = concentration of ribosomes in polysome; M t = concentration of ribosomes at termination codon site; P = concentration of protein; E = concentration of free start/initiation codon site(s); k A = rate constant for free ribosome occupying the start codon site; k I = rate constant for ribosome entering the polysome pool; k E = rate constant for departure of ribosome from polysome pool; k T = rate constant for release of synthesized protein (and release of ribosome); k d = degradation rate constant of protein product. (MR) = total concentration of mRNAs, a constant (not shown).

Fig. 3 .

 3 Fig. 3. Schematic of model 2 This figure depicts a schematic of model 2. This schematic is similar to that of model 1, except that the role of release factor protein P in regulating its own termination is shown. The variables are similar to those of model 1.

Fig. 4 .

 4 Fig. 4. Schematic of model 3 This figure depicts a schematic of model 3, indicating the presence of a premature stop codon. There are now two polysome pools, and two possible termination sites. At the premature stop codon, translation can either get terminated, or readthrough can occur, mediated by a tRNA species indicated by T. RP 1 = concentration of ribosomes in first polysome pool; M t1 = concentration of ribosomes at the premature termination codon site; P t = concentration of truncated protein; T = concentration of tRNA; (RP 2 ) = concentration of ribosomes in second polysome pool; P = concentration of protein; k E1 , k E2 = rate constants for exit from polysome pools 1 and 2 respectively; k t1 = (second order) rate constant for readthrough of the premature stop codon; k T a , k T b = rate constants for release of protein at the premature and final termination sites respectively respectively; k d = degradation rate constant of protein product based on experimentally measured half-life; k dt = degradation rate constant of truncated protein product.

Fig. 5 .

 5 Fig. 5. Schematic of model 4 This figure depicts a schematic of model 4, indicating the presence of a premature stop codon along with feedback in termination (the full length protein is involved in termination events at both termination sites). The variables are similar to that in model 3; note that in model 4, k T a and k T b are second order rate constants. As in model 3, readthrough of the premature stop codon can occur mediated by the tRNA.

Fig. 6 .

 6 Fig. 6. Simulation results from model 1. Temporal evolution of concentrations of (a) polysome pool (RP ) and (b) protein (P ) from model 1 simulations. Both eventually reach steady states. Parameter values: k A =400, k I =0.1, k E =0.0022, k T =0.1, k d =1.86e-4 all s -1 ; (MR)=1.26e-4 moles/volume; r 0 (total concentration of ribosomes) =1.26e-3 moles/volume; P 0 (initial concentration of protein) =1.6e-4 moles/volume.

Fig. 7 .

 7 Fig. 7. Simulation results from model 2. Temporal eveolution of concentrations of (a) polysome pool (RP ) and (b) protein (P ) from model 2 simulations. k A =400, k I =0.1, k E =0.0022, k d =1.86e-4 all s -1 ; k T =10/mole.s; (MR)=1.26e-4 moles/volume; r 0 =1.26e-3 moles/volume; P 0 =1.6e-4 moles/volume.

Fig. 8 .

 8 Fig. 8. Simulation results from model 3. Simulated concentrations of ribosomes in (a) ribosome pool one (RP 1 ),(c) ribosome pool two (RP 2 ), and that of (b) truncated protein (P ) and (d) full length protein (P t ) for model 3.k A =400, k I =0.1, k E1 =0.007, k t1 T =100, k E2 =0.0021, k d =1.86e-4, k dt =9.6e-5, k T a =0.1, k T b =0.1 all s -1 ; (MR)=1.26e-4 moles/volume; r 0 =1.26e-3 moles/volume; P 0 =1.6e-4 moles/volume. The premature stop codon is located 100 codons along the mRNA.

Fig. 9 .

 9 Fig. 9 .Effect of tRNA concentrations in model 3 . The effect of tRNA concentrations in determining readthrough and its effect on concentrations of full length and truncated proteins is studied here. Log-log plots of steady state concentrations for full-length protein, P (circles) and truncated protein, P t (crosses) versus a scaled product of the rate constant k t1 multiplied by the tRNA concentration in model 3 is shown. Results obtained by simulations are shown in (a) and analytical results are shown in (b).In general, we see that higher tRNA concentrations lead to higher concentrations of full length

Fig. 10 .

 10 Fig. 10. Simulation results from model 4. Typical concentrations of ribosomes in (a) ribosome pool 1 (RP 1 ), (c) ribosome pool 2 (RP 2 ), and of (b) truncated protein (P ) and (d) full length protein (P t ) for model 4 simulations. All common parameters are as in the previous figure; k T a =10 1/mole.s;

  what the roles of additional features such as feedback effects and the presence of premature stop codons are. What follows is a model-by-model description of these results and examination of the steady state analytical solutions.

1 Wolfram Research, Inc., Mathematica, Version 7.0, Champaign, IL (2008). 2 CRNT, Version 1.02 can be freely downloaded from Martin Feinberg's web page at www.che.eng.ohio-state.edu/feinberg/crnt A c c e p t e d m a n u s c r i p t

  α 1 .P ; (RP 2 ) = α 2 .P ; M t1 = α 3 .P ; (RP 1 ) = α 4 .P

		(48)
	P t = α 5 .P ; M o = α 6 .P ; r.E = α 7 .P	(49)

Acknowledgement

We gratefully acknowledge funding from EPSRC (co-funded with BBSRC), UK