Table Of Content

NonamemanuscriptNo. (willbeinsertedbytheeditor) Virus Replication as a Phenotypic Version of Polynucleotide Evolution DedicatedtothememoryofourdearcolaboratorandfriendFranciscodeAssisRibasBosco. FernandoAntoneli⋆ FranciscoBosco · · DiogoCastro LuizMarioJanini · 3 January11,2013 1 0 2 Abstract Inthispaperwerevisitandadapttoviralevolutionanapproachbasedonthethe- n oryofbranchingprocessadvancedbyDemetrius,SchusterandSigmund(“Polynucleotide a evolutionandbranchingprocesses”,Bull.Math.Biol.46(1985)239–262),intheirstudyof J polynucleotide evolution. Bytaking intoaccount beneficial effects weobtain anon-trivial 0 multivariategeneralizationoftheirsingle-typebranchingprocessmodel.Perturbativetech- 1 niquesallowsustoobtainanalyticalasymptoticexpressionsforthemainglobalparameters ] of themodel which lead tothe following rigorous results: (i) anew criterion for“nosure E extinction”,(ii)ageneralizationandproof, forthisparticularclassofmodels,ofthelethal P mutagenesiscriterionproposedbyBull,SanjuańandWilke(“Theoryoflethalmutagenesis . o forviruses”,J.Virology18(2007)2930–2939), (iii)anewproposalforthenotionofrelax- i ationtimewithaquantitativeprescriptionforitsevaluation,(iv)thequantitativedescription b of the evolution of the expected values in in four distinct “stages”: extinction threshold, - q lethalmutagenesis,stationary“equilibrium”andtransient.Finally,basedonthesequantita- [ tiveresultsweareabletodrawsomequalitativeconclusions. 3 Keywords Viralevolution Branchingprocesses Phenotypicmodel v · · 3 MathematicsSubjectClassification(2010) 92Dxx 60J80 60J85 5 · · 3 6 ⋆correspondingauthor. . 4 FernandoAntoneli 0 DepartamentodeInforma´ticaemSau´de,UniversidadeFederaldeSaõPaulo,SaõPaulo,SP,Brazil. 2 E-mail:[email protected] 1 FranciscoBosco : Laborato´riodeBiocomplexidadeeGenoˆmicaEvolutivaandDepartamentodeMicrobiologia,Imunologiae v Parasitologia,UniversidadeFederaldeSaõPaulo,SaõPaulo,SP,Brazil. i X E-mail:[email protected] r DiogoCastro a DepartamentodeMedicina,UniversidadeFederaldeSaõPaulo,SaõPaulo,SP,Brazil. E-mail:[email protected] LuizMarioJanini DepartamentodeMicrobiologia, ImunologiaeParasitologiaandDepartamentodeMedicina,Universidade FederaldeSaõPaulo,SaõPaulo,SP,Brazil. E-mail:[email protected] 2 Antoneli,Bosco,CastroandJanini 1 Introduction RNAvirusesexhibitapronouncedgeneticdiversity(Domingoandetal.,1985).Thisvari- abilityallowsRNAvirustobetteradapttoenvironmentalchallengesasrepresentedbyhost and therapy pressures (Domingoetal, 1998). Due to the lack of a proofreading activity of viral RNA polymerases (average error incorporation rate of order 10 4 per nucleotide, − perreplication cycle(Steinhaueretal,1992)),short generation timesand hugepopulation numbers,RNAviralpopulationsmaybeviewedasacollectionofparticlesbearingmutant genomes. As a consequence of high mutation rates, frequencies of mutants depend not only on their level of adaptation but on the probability of being faithfully replicated during viral genomic RNA synthesis. Several studies (Domingoetal, 1998; ElenaandSanjuań, 2005) havelookedatviraldiversificationprocessesasacontributingcauseofdiseaseprogression andoftherapeuticstrategiesshortcomingsincludingvaccinetrials.Ithasbecomeimportant tounderstandtheprocessbywhichvirusacquirediversityandthedynamicsandfluctuations ofthisdiversityintime.However,understanding viralevolutioninvivohasproventobea veryunwieldyaccomplishmentduetothehugenumberofvariablesinvolvedintheinterplay betweenvirusesandtheirhosts. Forthelast30years,quasi-speciestheoryhasprovidedapopulation-based framework tomodelRNAviralevolution.Quasi-speciestheoryisamathematicalframeworkthatwas initially formulated to explain the evolution of macromolecules (Eigen, 1971). More re- cently, it has been used to describe the evolutionary dynamics of RNA viruses. However, quasi-species theory is based on deterministic equations and so has a number of serious drawbackswhenappliedtorealisticexperimentalsystemssincethereareimportantsources ofstochasticitythatmustbetakenintoaccount. Traditionally, in an effort to make the viral evolution process more palpable, several groups have addressed this subject from different points of view. There is a substantial amount ofpublications thatstudiedviruspopulations during theirevolution inexperimen- talsettings,forinstance,cellcultures(Escarm´ısetal, 2006),bychallengingtheviruswith population bottlenecks (Ojosnegrosetal, 2010b,a), or the introduction of antiviral drugs (Crottyetal, 2001), including mutagens, or another competing viral population. Several othergroups havestudiedtheprocessofviral evolution awayfromthebench, usingcom- putational and mathematical tools (Wilkeetal, 2001; Bulletal, 2007). The main lesson learnedfromtheseeffortsisthatinordertoescapefromoversimplifying theinterplaybe- tween virus and hosts, one needs to take into account a few hard rules based on the ex- perimentaldatathathasbeengeneratedbythecommunityofinvestigatorsaddressingviral evolution. In this paper we revisit and adapt to viral evolution an approach based on the theory of branching process advanced by Demetriusetal (1985) in their study of polynucleotide evolution. In fact, almost any stochastic model of asexual replication has an underlying branchingprocess,whichremainsimplicitandundefinedinmostofthestudiesinthissub- ject.Forinstance,(Manrubiaetal,2003;AguirreandManrubia,2008;Aguirreetal,2009; Cuestaetal,2011;Capitańetal,2011;Cuesta,2011)haveinvestigatedprobabilisticmodels in the form of a linear mean field approximation without explicitly defining the underly- ingstochasticprocessmodelingthemicroscopicdynamicsofparticlereplication.However, thereareafewofthem(Demetriusetal,1985;Hermissonetal,2002;Wilke,2003)thattake adifferentpathandexplicitlydefinethestochasticprocessinordertobringthemathemat- icaltheoryofbranchingprocessestobear.Thisattitudehassomevirtues,sinceitprovides VirusReplicationasaPhenotypicVersionofPolynucleotideEvolution 3 powerfultools,thathavebeenperfectedinthepastseveraldecades,allowingonetoextract quantitativeresultsonarigorousbasisinclearconceptualframework. Inthepresentworkwepursuethesamepathas(Demetriusetal,1985;Hermissonetal, 2002;Wilke,2003),byputtingmoreemphasisonthestudyofanexplicitgeneratingfunction oftheunderlyingbranchingprocessandapplyingthemathematicaltoolsasameanstoget newinsightsandperspectivesonthedynamicsofevolvingviruspopulations.Beforesetting upourmathematicalframeworkletusgiveabroadoverviewofthescenerythatweconcieve asthecontextofthemodelsdescribedhere. Viral infections may result in different clinical outcomes which are consequences of interactionsbetweenthevirusanditshost.Theconstantarmsracebetweenhostsandvirus caneithereliminateorperpetuateviralinfections.Someviraldiseasesareselflimitedending withtheeradicationoftheinfectingviruswhileothersleadtowardspersistentinfections. Duringitsadaptationtonewhostsviruspopulationswillshowanintensefluctuationin sizeandstructureinresponsetotheselectivepressuresimposedbythehostdefenses.Soon aftertheinfectionvirusparticlesstarttocolonizethenewenvironmentandwillseekwithin thehost forthecelltypes that willbettersupport theproduction ofanexuberant progeny. Atthismomentandwiththetotalabsenceofanimmuneresponsetheviralpopulationmay experienceanexponential growth.Thisgrowthisclinicallyreflectedbytheappearanceof acutesymptoms and theemergence ofa peak of viremiadetected inthe peripheral blood. Theemergenceofselectivepressuresrepresentedbythecombinationofcomponentspresent intheinnateandadaptiveresponsesoftheimmunesystemwillcauseadropintheviralload tolevelswellbelowtheinitialviremiavalue.Atthispointthevirusandhostdefensesmay reachastationary“equilibrium”referredtoastheviralsetpointinsomechronicinfections asthe onecaused bythe human immunodeficiency virus (HIV). On theotherhand, ifthe burdenimposedbythehostpressuresovercomesthevirussurvivalcapacity,theviralpopu- lationwillbecomeextinct.Theinitialexponentialgrowthandsubsequentpopulationsizere- ductionexperiencedbytheincomingviruscanbereferredtoasthetransientphaseofaviral growthcurve.Thetransientphasewillrepresenttheviruspopulationfluctuationsduringthe periodfromitsintroductioninthenewhostuntilthemoment theviralsetpointisreached. Inpersistentinfectionswithoutmajoroscillations,inthehostselectivepressuresorinvirus adaptability, virus populations will show an asymptotic behavior and will perpetuate in a stationaryequilibriumregime.InviruswithgenomemoleculesrepresentedbyRNAwhich are able to tolerate high mutational rates, the stationary equilibrium is actually a balance betweenthetwomajorforcesinvirusevolution:mutationandselection.Asamatteroffact, thetransient phaseofavirus growth hasbeen alsocalledthe recoverytimeasanallusion to virus populations recovering from bottleneck passages that correspond todrastic population sizereductions. When selectivepressures against the incoming virus aretoo strong forthevirustocouplewiththemthestationaryphaseofthevirusgrowthisneverattained andtheasymptoticbehaviorisabsent.However,changesthatwillcausearecrudescencein thehostselectivepressuresmaydislocatevirusfromthestationaryequilibriumtowardsits extinction. Therapeutic strategies as antiviral drugs and therapeutic vaccines may precipi- tateviruspopulationsinregionsofinstabilityrepresentedbyhostileadaptationsmarkedby subsequentdropsinvirusreplicationcapacityandprogenysizes.Thisphaseofviralgrowth namedhereasextinctionthreshold representsamoment ofuncertaintywhenthevirusex- tinctionisstronglypossiblebutunpredictable.Movingforward,whenviruspopulationsare unabletocounteractthedeleteriouseffectsofallnaturalandchemicalpressurescombined, andaveragevirusprogenydropstoavaluebellowoneviableparticle,theprocesscalledthe lethalmutagenesisisstarted.Atthismomenttheviruswillbecomeextinctalmostcertainly. Virusextinctionscanbeseenduringselflimitedinfectionsevenbeforeanasymptoticbehav- 4 Antoneli,Bosco,CastroandJanini iorisreachedorwhenchronicinfectionsareclearedduetooptimizationofhostpressures ortotheadditionofantiviralchemicals. Once the mathematical framework is established, the scenary described above raises severalquestionsabouttheevolutionofaviralpopulation withinitshostcanbeposedand answeredinapreciseway.Inthisworkweshalladdressthefollowingquestions: (i) Whatistheaveragenumberofparticlesinthepopulationineachgeneration? (ii) Whatistheprobabilityofextinctionofthepopulationwithinitshost? (iii) Howdoesthe“transition”toextinctionhappens? (iv) Whenthepopulationdoesnotbecomeextinct,whatistheasymptoticbehaviour? (v) Howlongittakesforthepopulationtoreachitsasymptoticstate? In order to answerthese questions we propose anon-trivial multivariate generalization of thesingle-typebranchingprocessofDemetriusetal(1985). Byemploying perturbative techniques weareabletoobtain analytical expressions for themainglobalparameterofthemodel,namelytheaveragegrowthrate,representedhereby themalthusianparameter ofthebranchingprocess.Ourresultscorroboratewiththestudy ofBulletal(2007)byshowingthatthesufficientconditionforlethalmutagenesisinvolves mutational andecological aspects.Theyarrivedataconjectural criteriaforlethalmutage- nesisbyaheuristicandintuitiveapproachofgreatgeneral applicability.Fortheparticular classofmodelsconsideredhereitispossibletostateandrigorouslyproveanextinctioncri- terionwhichis,underaproperinterpertation,equivalentetothelethalmutagenesiscriterion ofBulletal(2007).Inotherwordsthebranchingprocessformalismcanprovideanewand interesting perspective on this problem. Weobtain anew criterionfor non-extinction of a viralpopulationanddescribefourdistinct“stages”ofevolutionofaRNAviruspopulations: extinctionthreshold,lethalmutagenesis,stationary“equilibrium”anddecayofthetemporal auto-correlation. 2 PhenotypicModelsforViralEvolution In this section we describe a model for viral evolution that is naturally represented by a multivariatebranchingstochasticprocessgeneralizing, inanon-trivialway,thesingle-type branchingprocessstudiedbyDemetriusetal(1985).Forthesakeofmotivationwestartby recallingaprobabilisticmodelintroducedbyLa´zaroetal(2002). Weinterpret thenotion ofmutation probability as ageneral effect of probabilistic na- turewithdirectimpactonthereplicationcapabilityofindividualviralparticles,considered hereasameasureoftheparticle’sfitnesscharacterizing itsphenotype. Thiseffectissum- marized by the definition of astationary probability distribution whichis used toset up a Galton-Watsonbranchingprocess(WatsonandGalton,1874)forthetemporalevolutionof the viral population. This probability distribution gives appropriate parameters to classify theasymptoticbehavioroftheviralpopulationandtodescribesomeofthenon-equilibrium propertiesofthemodel. In other publications on the same subject the concept of mutation is extensively used as the cause of replication capacity change. Understanding that those changes constitute anobservableoutputduetomanydifferentfactors(ofgeneticandnon-geneticnature),we prefertousethegeneralterm“effect”overthereplication capacitytocharacterizethethree possiblechanges(deleterious,beneficialandneutral)thatmayhappenwiththeviralparticle whenitreplicates. VirusReplicationasaPhenotypicVersionofPolynucleotideEvolution 5 2.1 TheProbabilisticModel A number of viral infections starts with the transmission of a relatively small number of viralparticlesfromoneorganismtoanotherone.Theinitialviralpopulationstartsreplicat- ingconstrainedbytheunavoidable interactionwiththehostorganismandevolves intime towards aneventual equilibrium. Eachparticle composing thepopulation replicates inthe cellularcontextthatmaydifferfromcelltocell.Moreovereachparticlehasdifferentreplica- tioncapabilitiesduetothenaturalgenomicdiversityfoundinviralpopulationsingeneral. Therefore, it is reasonable to consider the viral population as a set of particles divided in groupsofdifferentreplicationcapabilitiesmeasuredintermsofthenumberofparticlesthat oneparticlecanproduce.Themodelsconsideredheredonottakeintoaccountanyinforma- tionaboutthegenomicdiversityofanyreplicatingclassandthereforeitshouldbeclassified asphenotypicmodel. Considerthatthewholesetofparticlescomposingtheviralpopulationreplicatesatthe sametime in such away that theevolution of the population is described as a succession of discrete viral generations. This assumption crucially depends on the cleardefinition of the time needed for a particle to replicate, referred by virologists as generation time. As it depends on the cell environment it is clearthat this time period may vary from particle toparticle, replicating indifferent cells,insuch away that ameaningful concept isadis- tribution of replication times with a possible well defined mean value. The dispersion of thereplicationtimescanbeconsidered smallifwerestrict ourselves tohomogeneous cell populations. Undertheseconditions, onemay considerthat noparticlecan bepart oftwo successivegenerations,thatis,thegenerationsarediscreteandnon-overlapping. Supposethatapopulation ofvirusesthatstartevolving fromaninitialsetofparticles, which is partitioned into classes according to their mean replication capacity, that is, for each integer r=0,...,R there is a class of particles labeled by r and a random variable assumingnon-negativeintegervalueswhoseprobabilitydistributiont (k)satisfiesE(t )=r r r andt (k)=d . 0 k0 Inasmuch as the process of replication is controlled by chemical reactions involving specificenzymes and thetemplate, itisreasonable toassume ameanbounded replication capacityperparticlethatispossiblytypicalforeachspecificvirus.Hencethereisamaxi- mummeanreplicationcapacityRimposedbythenatural limitingconditions underwhich anyparticleofthepopulationreplicates. In the process of replication of a viral particle errors may occur at each replication cycle in the form of point mutations with possible impact on the replication capacity of theprogeny particles.Duetotheintrinsicstochasticcomponent ofchemicalreactionsitis naturaltotreatthispointmutationalcauseasprobabilistic.Anotherpossiblecauseofchange inthereplicationcapabilityintheviraloffspringisclearlyrelatedtothecellularenvironment wherethereplicationprocesstakesplace.Asaresultthetimeevolutionofviralpopulations shouldbeviewedasaphysicalprocessstronglyinfluencedbystochasticity.Therefore,the combined action of genetic and non-genetic causes may produce basically three types of replicative effects with associated probabilities, at the particles scale, applicable to every singlereplicationevent: – deleterious effect d: the mean replication capacity of the copied particle decreases by one.Notethatwhentheparticlehascapacityofreplicationequalto0itwillnotproduce anycopyofitselfontheaverage. 6 Antoneli,Bosco,CastroandJanini – beneficial effect b: the replication capacity of the copy increases by one. If the mean replicationcapacityisalreadythemaximumallowedthenthemeanreplicationcapacity ofthecopieswillstaythesame. – neutraleffectc:themeanreplicationcapacityofthecopiesremainthesameasthemean replicationcapacityoftheparentalparticle. Theonlyconstraintthesenumbersshouldsatisfyisb+c+d=1.Inthecaseofinvitroex- perimentswithhomogeneouscellpopulationstheparametersc,dandbmaybeconsidered essentiallyasmutationprobabilities. The assumption that b is very small when compared with d and c is justified by several experimental results.Thefrequencies betweenbeneficial, deleterious andneutral mutations appearing in areplicating population have been already measured by prior studies (Mirallesetal, 1999; ImhofandSchlotterer, 2001; KeightleyandLynch, 2003; Orr, 2003; Sanjuańetal, 2004; Carrascoetal, 2007; Eyre-WalkerandKeightley, 2007; Pareraetal, 2007; Rokytaetal, 2008). Taking their results together, it is reasonable to conclude that beneficial mutations could be as low as 1000 less frequent than eitherneutral or deleterious mutations. As a result, the viral population would be submitted to a large number of successivedeleteriousandneutralchangesandacomparativelysmallnumberofbeneficial changes.Theparticularcasewheretherearenobeneficialeffectsintime(b=0)isinterest- ingnotonlybecauseofbiologicalreasonsbutalsoduetoamathematicalproperty,namely, thespectralproblemassociatedtoitsmeanmatrixis“completelysolvable”.Thisproperty will be crucial for us to implement the perturbative expansions in the general case where b=0. 6 2.2 TheMultitypeBranchingProcess Instead of working directly with a probabilistic model we will introduce a new generat- ingfunction,whichisanon-trivialmultivariategeneralizationofthesingle-typebranching processproposedbyDemetriusetal(1985).Forfullpresentationsofthetheoryofbranch- ingprocessessee(Harris,1963;AthreyaandNey,1972;KimmelandAxelrod,2002).Brief accounts of the theory branching process, highlighting the main results weneed here, are providedby(Demetriusetal,1985,Sec.2)andAntonelietal(2013a). Inordertofindourprobabilitygeneratingfunctionweobservethatthereplicationpro- cessissimplyaBernoullitrialwiththreepossibleoutcomes:anewlyproducedparticlemay haveenduredadeleterious,neutralorbeneficialeffect.Therefore,theoffspringdistribution shouldbeatrinomialdistribution(seeFeller(1968))andtheprobabilitygeneratingfunction ofthephenotypicmodelisf =(f ,f ,...,f )withcomponents 0 1 R ¥ f (z ,z ,...,z )= (cid:229) t (k)(dz +cz +bz )k. (1) s 0 1 R s s 1 s s+1 − k=0 with“consistencyconditions” f (z ,z ,...,z )=1andz =z . 0 0 1 R R+1 R Noteitistrivialtofurthergeneralizethemodelinordertoinclude“higherorderreplica- tive effects” that changes the mean replication capacity by 2, 3,... and thus replacing ± ± thetrinomialdistributionbyamultinomialdistribution.Itisalsoworthtomentionthatthere areotherpossiblevariationsofthesemodelsthatsharethesameessentialpropertiesandare moreadequateindifferentcontexts: VirusReplicationasaPhenotypicVersionofPolynucleotideEvolution 7 WithZeroClass: Inthisversionthemodelnaturalyfitsthesymmetryofthebinomialdis- tributionsinceitcontainstheparticlesofclassr=0whicharegeneratedbytheparticles fromclassr=1. WithoutZeroClass: Inthisvariation,theparticleclass0isomittedandthustheprobabil- itygenerating function has Rvariables and Rcomponents: setthevariable z =1and 0 omitthefirstcomponent f .Particlesofclassr=1undergoingadeleteriouschangeare 0 eliminatedinthenextgeneration.Inthisformulationthemodelispositivelyregular. Withtheconventiondescribedabove,itiseasytoseethatintheone-dimensional case withb=0oneobtainsthefollowinggeneratingfunction ¥ ¥ f(z) = (cid:229) t(k)((1 c)+cz)k = (cid:229) t(k)(1 c(1 z))k. (2) − − − k=0 k=0 Thisisexactlythegeneratingfunctionofthesingle-typemodelproposedby(Demetriusetal, n 1985,p.255,eq.(49))fortheevolutionofpolynucleotides.Intheirformulation, c=p is theprobabilitythatagivencopyofapolynucleotideisexact,wherethepolymerhaschain lengthofn nucleotidesandthereisafixedprobability pofcopyingasinglenucleotidecor- rectly.Thereplicationdistribution t providesthenumberofcopiesapolynucleotideyields beforeitisdegradedbyhydrolyses. 2.3 EvolutionoftheMeanValues We shall introduce the notation Zr =1 for the condition Z =e , which is the initial 0 0 r population consisting of one particle of class r and no particles of other classes. Thus P(Z =iZr =1)=z (i).Abasicassumptioninthetheoryofbranchingprocessesisthat 1 | 0 r all the first moments are finite and that they are not all zero. Then one may consider the meanevolutionmatrixorthematrixoffirstmomentsM = M whichdescribeshowthe ij { } averages Z ofthesub-populationsofparticlesineachreplicationclassevolvesintime: n h i Z =MnZ or Z =M Z . (3) n 0 n n 1 h i h i h − i ThemeanmatrixM = M isdefinedas ij { } ¶ f M =E(Zi Zj=1) or M = j(1), ij 1| 0 ij ¶ zi fori,j=0,...,Rand1=(1,1,...,1). Inordertocompute themeanmatrixof(1)weobservethat only thepartial derivative of a component f with respect to z , z and z are non-trivial, the others are zero. r r 1 r r+1 Therefore,themeanmatrixis − 0 d 0 0 0 0 ··· 0 c 2d 0 0 0  ···  0 b 2c 3d 0 0 M = ··· (4) 0 0 2b 3c 0 0  ... ... ... ... ·.·.·. (R 1)c Rd   −  0 0 0 0 (R 1)b R(c+b)  ··· −    Thusshowingthatourgeneratingfunctionindeedprovidesanon-trivial(i.e.,non-trianglar) generalizationof(2).Interestingly,themeanmatrixM providedbythisclassofmodelsis 8 Antoneli,Bosco,CastroandJanini atridiagonalmatrix,whichisanubiquitoustypeofmatrixappearinginseveralfieldsrang- ingfromstatisticalsignalprocessingGray(2006),informationtheoryGrenanderandSzego¨ (1958),latticedynamicalsystems(Antonelietal,2005). Themeanmatrixobtainedabovecoincidewithsomeofthematricesinthelinearmean field approximation studied in Manrubiaetal (2003); AguirreandManrubia (2008), Aguirreetal(2009);Cuestaetal(2011);Capitańetal(2011);Cuesta(2011).However,itis importanttostressthattheseworksconsiderthatthetotalpopulationZ evolvesbymulti- n plicationbythemeanmatrix,whilehereitistheevolution oftheaveragevalue Z that n h i isdescribedbythemeanmatrix.However,themeanmatrix M dependsontheprobability distributionst (k)onlythroughitsmeanvaluer.Therefore,thereareinfinitelymanydistinct r branching processes withthe samemeanmatrix and theonly way totellthem apart isby lookingattheirfluctuationsorsecondmomentproperties. Letr (M)denotethespectralradiusofthemeanmatrixM,thatis,ifl 6 6l are 1 R theeigenvaluesofMthenr (M)=l .SinceMisanon-negativematrix,ithas·a·t·leastone R largest non-negative eigenvalue which coincides with its spectral radius (see Gantmacher (2005) or Meyer (2000)). When the largest eigenvalue is positive we shall call it, follow- ingKimmelandAxelrod(2002),themalthusianparameterm=r (M)=l ofthebranch- R ingprocess(seealsoJagersetal(2007)). Themainclassificationresultofmultitypebranchingprocesses,inthepositivelyregular case,statesthatthereareonlythreepossibleregimes(seeHarris(1963)orAthreyaandNey (1972)): sub-critical (m<1), super-critical (m>1) and critical (m=1). Note however, that when b=0, the corresponding branching process is not positively regular. Neverthe- less, there is a generalization of the classification of multitype branching processes, due toSevastyanov(seeHarris(1963); Jiˇrina(1970)), thatallowsus toinclude thecaseb=0 withinthesamethreeregimesasabove. 3 MalthusianParameterofthePhenotypicModel Inthissectionweshallemploytheperturbationtheoryforeigenvaluesandeigenvectorsin ordertostudytheaveragedynamicsofthephenotypicmodel,withbsmall.Inordertomain- tainthemainflowofthetextwehaveomittedmostofthecomputationsandestimateswhich have been placed in the appendix for the convenience of the reader. See also Wilkinson (1965)orThompsonandMcBride(1974)formoredetails. Methods of spectral perturbation have been used in the study of evolution of macromolecules(SwetinaandSchuster,1982;ThompsonandMcBride,1974),wherethethestart- ing point is a diagonal matrix and the perturbation adds other off-diagonal elements. Our approachissimilar,butourstartingpointisanuppertriangularmatrixandthefinalmatrix istridiagonal. In order to implement this program it is necessary to know all the eigenvalues of the unperturbed matrix and their associated left and right eigenvectors. The idea is to write themeanmatrix(4)assumofamatrix whosespectral problem isexactlysolvablewitha “perturbation matrix”depending onasmallparameter, insuchawaythattheperturbation vanishesastheparametergoestozero. LetM bethemeanmatrixof(4)thephenotypic model.Then,sincec=1 b d we − − maywrite M =M +bP, 0 whereM isthemeanmatrixofthesimplephenotypimodel: 0 VirusReplicationasaPhenotypicVersionofPolynucleotideEvolution 9 0 d 0 ... 0 0(1 d) 2d ... 0  −  0 0 2(1 d) ... 0 M0=0 0 0− ... 0  (5)   ... ... ... ... Rd    0 0 0 ... R(1 d)  −    andPisthematrix 0 0 0 0 0 ... 0 0 0 1 0 0 0 ... 0 0  −  0 1 2 0 0 ... 0 0 − 0 0 2 3 0 ... 0 0 P =0 0 0 −3 4 ... 0 0. (6)  −  ... ... ... ... ... ... ... ...   0 0 0 0 0 ... (R 1) 0  − −  0 0 0 0 0 ... (R 1) 0  −    It iseasytoseethatM isthemean matrix ofthephenotypic model when onesets b=0 andc=1 d inthegeneratingfunction(1).Aninterestingfeatureofthisparticularcaseis − thatitsspectralproblemiscompletelysolvableasfunctionsoftheparametersdandR.The eigenvaluesl 0ofthemeanmatrixM are: r 0 l 0=rc=r(1 d) r=0,...,R. r − Inparticular,themalthusianparameteristhelargestpositiveeigenvalue m0=l 0=Rc=R(1 d). (7) R − Also,notethattheoperatornormofP is P =2(R 1). Nowwewritetheeigenvaluel ofMkaskafuncti−onoftheparameterb,expandedasa r powerseries l =l 0+l 1b+l 2b2+..., r r r r wherel 0isthecorrespondingeigenvalueofM .Clearly,l (b) l 0asb 0.Thehigher order perrturbation terms l i, (i=1,2,3,...) are0written in rterms→ofrall the→eigenvalues l 0 r s (s=0,...,R)ofM andtheirassociatedleftandrightnormalizedeigenvectorsv0andu0. 0 s s WeareinterestedinthemalthusianparametermofthematrixM: m=m0+m1b+m2b2+..., wherem0=(1 d)R.Alengthycalculationgivesthesecondorderexpansionofthemalthu- − sianparameterm: (R 1)d m=(1 d)R+ − Rb − 1 d − (8) (1+dR/2)(R 1)2d − Rb2+O(b3). − (1 d)3 − Theparameterspaceofthemodelis (b,d,R) R2 N b+d61 .Ifweforgetthe { ∈ × | } discreteparameterR,thenthissetcanbeidentifiedwiththetriangle (0,0),(0,1),(1,0) in { } the(d,b)-plane(seeFIG.1). 10 Antoneli,Bosco,CastroandJanini Fig.1 Parameterspace(d,b)ofthephenotypicmodelwiththecriticalcurvesm=1forR=2,3,4. Nowwecanmakesomequantitativepredictionsaboutthebehaviourofthemodel. One can localize in the parameter space (b,d) each of the three regimes of the phenotypic model (sub-critical, super-critical, critical) by considering, for each fixed R, the criticalcurvem(b,d,R)=1(seeFIG.1).Theregionofsub-criticality(foreachfixedR)lies on the right sideof thecritical curve and the region of super-criticality (foreach fixed R) liesontheleftsideofthecriticalcurve. Anotherinterestingfeatureoftheparameterspaceisrelatedtothepointsofintersection ofthecriticalcurvem(b,d,R)=1withtheboundarycurvesb+d=1andb=0.Thepoints ofintersectionwiththeboundary curve b=0aregiveby d (R)=1 1/Randthepoints c − of intersection with the boundary curve b+d =1 have non-zero coordinates and can be writtenas(1 b ,b )whereb (R)istheb-coordinateand1 b (R)isthed-coordinate.For c c c c instance,itis−easytofindtheexact valueofbc forR=2:−bc(2)=1−√12 ∼=0.2928. This calculationgivesanimportantprediction:ifb>0.2929thenthepopulationcannotbecome extinctbyincreasingdtowardsitsmaximumvalue(=0.7070),withprobabilitygoingto1 ∼ asRbecomeslarge. Moregenerally,onehasthefollowingcriterionofnosureextinction.Considerthephe- notypicmodelwithparameters (b,d,R),withRlargeenough.Ifb=0issufficientlysmall and(uptoorderO(b2)) 6 bR(R 1)2>1 (9) − then,asymptoticallyalmostsurely,thepopulationcannotbecomeextinctbyincreasingthe deleteriousprobabilitytowardsitsmaximumvalue.Inotherwords,whentheinequality(9) issatisfiedtheprocess canhave onlyoneregime,namely, the supercritical, whichrenders themodelcompletelytrivial,sincethereisnotransition.