The Abridgment and Relaxation Time for a Linear Multi-Scale Model Based on Multiple Site ...

19 Pages·2015·0.82 MB·English

RESEARCHARTICLE
The Abridgment and Relaxation Time for a
Linear Multi-Scale Model Based on Multiple
Site Phosphorylation
ShuoWang,YangCao*
DepartmentofComputerScience,VirginiaTech,Blacksburg,VA,UnitedStatesofAmerica
* [email protected]
Abstract
Randomeffectincellularsystemsisanimportanttopicinsystemsbiologyandoftensimu-
a11111 latedwithGillespie’sstochasticsimulationalgorithm(SSA).Abridgmentreferstomodel
reductionthatapproximatesagroupofreactionsbyasmallergroupwithfewerspeciesand
reactions.Thispaperpresentsatheoreticalanalysis,basedoncomparisonofthefirstexit
time,fortheabridgmentonalinearchainreactionmodelmotivatedbysystemswithmultiple
phosphorylationsites.Theanalysisshowsthatiftherelaxationtimeofthefastsubsystemis
muchsmallerthanthemeanfiringtimeoftheslowreactions,theabridgmentcanbeapplied
OPENACCESS withlittleerror.Thisanalysisisfurtherverifiedwithnumericalexperimentsformodelsofbis-
tableswitchandoscillationsinwhichlinearchainsystemplaysacriticalrole.
Citation:WangS,CaoY(2015)TheAbridgmentand
RelaxationTimeforaLinearMulti-ScaleModel
BasedonMultipleSitePhosphorylation.PLoSONE
10(8):e0133295.doi:10.1371/journal.pone.0133295
Editor:RamonGrima,UniversityofEdinburgh,
UNITEDKINGDOM Introduction
Received:February27,2015 Withtherapiddevelopmentinsystemsbiology,biologicalmodelshavebecomemoreand
Accepted:June25,2015 morecomplex.Inmanycomplexbiochemicalsystems,stochasticeffectsincellsareofparticu-
larconcernbecausesomespecies,suchasgenesandRNAs,presentwithsmallcopynumbers
Published:August11,2015
inthosesystems.Inordertostudythestochasticeffects,stochasticbiochemicalmodelshave
Copyright:©2015Wang,Cao.Thisisanopen beenbuiltandsimulated.Gillespie’sstochasticsimulationalgorithm(SSA)[1,2]isoneofthe
accessarticledistributedunderthetermsofthe
mostimportantstochasticmethods.However,thecomputationalcostoftheSSAcanbevery
CreativeCommonsAttributionLicense,whichpermits
high,particularlyforsystemswiththemulti-scalefeature,whichhighlightsthescaledifferences
unrestricteduse,distribution,andreproductioninany
medium,providedtheoriginalauthorandsourceare amongreactions:Somereactionsfiremuchfasterthanothers,andthosefastreactionsoften
credited. quicklyreachequilibrium.SincetheSSAtrackseveryreactionfiring,inanSSAsimulation
mostcomputationalcostsarespentonfastreactions.Ontheotherhand,slowreactionsmay
DataAvailabilityStatement:Allrelevantdataare
withinthepaperanditsSupportingInformationfiles. bemoreimportantbecauseveryoftentheydrivethedynamicsofasystemwhenfastreactions
areinequilibrium.
Funding:YCaoreceivedfundingfromNational
Severalapproximationmethodshavebeenproposedtoutilizethismultiscalefeature.The
ScienceFoundationunderawardDMS-1225160and
CCF-0953590andtheNationalInstitutesofHealth keyideaoftheseapproximationmethodsistoapproximatethepropensitiesofslowreactions
underawardGM078989. bytakingaccountoftheeffectoffastreactionswhileavoidingfullsimulationofthem.Acom-
monpracticeistoreduceacomplexsystemtoasimplermodelthatwellapproximatestheorig-
CompetingInterests:Theauthorshavedeclared
thatnocompetinginterestsexist. inalsystem.Awell-knownexampleofsuchareductionistheMichaelis-Mentenequationof
PLOSONE|DOI:10.1371/journal.pone.0133295 August11,2015 1/19
TheAbridgmentandRelaxationTimeforaLinearMulti-ScaleModel
theenzyme-substratereactions,whichhasplayedanimportantroleinbiochemistry.Following
Gillespieetal.[3],wewillcallsuchareductionabridgment.Anabridgmentreplacesareaction
networkwithanetworkwithfewerreactionsandchemicalspecies.Anobviousadvantageofan
abridgmentisthereductioninthenumbersofreactionsandspeciesthatwehavetodealwith.
Anotheradvantagemightbespeedingupthenumericalsimulation.Butwealsohavetorealize
thedrawback.Anabridgmentisstillanapproximationbynature.Wehavetostudyconditions
thatanabridgmentisvalid.InGillespieet.al.[3],adetaileddiscussionwasgivenonthecondi-
tionthatthesystem
S Ðc1S !c3 S ð1Þ
1 2 3
c2
canbereducedtoasimpleone
S1!c S3; ð2Þ
wherecisgivensome“suitable”valuedependingonc ,c ,andc .Itwasconcludedthatthe
1 2 3
abridgmentisvalidifandonlyifoneofthefollowingconditionsissatisfied:
c (cid:2)c ; c (cid:2)c ; c (cid:2)c ; or c (cid:2)c : ð3Þ
2 1 3 1 1 3 2 3
Inthispaperweconsideramoregeneralproblem:Givenalinearchainreactionmodel
S Ðf1 S Ðf2 (cid:3)(cid:3)(cid:3)fÐn(cid:4)1S !fn S ; ð4Þ
1 2 n nþ1
b1 b2 bn(cid:4)1
whererateconstantsf’sandb’smaydependonothervariablesinthesystembutnotonS’s,
i i i
whatistheconditionthatthissystemcanbereducedto
S1(cid:4)!c Snþ1; ð5Þ
wherecdependsonthevaluesoff’sandb’s?
i i
System(4)isnotjustasimplegeneralizationof(1).Itrepresentsanimportantmodulein
biologicalsystemsaswell.Forexample,oneofthemostinterestingandwellstudiedbiological
systemsisthecellcyclemodel,whichaimstoaccuratelymodeltherepetitivesequenceof
eventsthatacellgrows,replicatesitscomponents,anddividesintotwodaughtercells.Akey
componentinthecellcyclemodelisbistableswitch,suchastheoneshowninFig1,where
Cdh1andClb2areregulativeproteinsinhibitingeachothertoformapositivefeedbackloop.
Onetypicalsettingofthetwoenzyme-substratereactionscanbetheMichaelis-Mentenrate
lawequation,whichassumesthattheconcentrationoftheenzymespeciesbemuchsmaller
thanthatofthesubstrate.However,thisassumptionmaynotalwaysbetrueincellcyclemod-
els.Particularlyinthisexample,sinceCdh1andClb2serveasenzymetoeachother.Thereis
nowaythattheaforementionedassumptionwillbevalidforbothdirections.Thusdifferent
kineticlawshavetobeadoptedtosolvethisdilemma.Elementaryreactions[1,2]areparticu-
larlydesiredinstochasticmodelingandsimulations.
Oneefforttosolvethisdilemmaistobreakdownsimpleswitchsystemsconstructedwith
Michaelis-Mentenratelawsintoelementaryreactionratelaws.Karetal.[4]managedto
replaceasimpleversionofthebuddingyeastcellcyclemodelwithaquitecomplexnetworkof
elementarychemicalreactionswithcorrespondingratelaws.Anotherinterestingdiscoveryby
Quetal.[5]showsthat,modelingtheregulativeproteinwithmultiplesitephosphorylation
mayalsoproducethe“nonlinearity”inthesystemandgeneratebistableswitchbehavior.The
ideaistoassumethattheregulativeproteinhavemultiplephosphorylationlevels,whichenable
PLOSONE|DOI:10.1371/journal.pone.0133295 August11,2015 2/19
TheAbridgmentandRelaxationTimeforaLinearMulti-ScaleModel
Fig1.Anexampleofbistableswitchinthecellcyclemodel.Eachdashedlineindicatestheenzyme
speciesononeenzyme-substratereaction.
doi:10.1371/journal.pone.0133295.g001
ittoundergomultiple(de)phosphorylationreactionsbetweendifferentlevels.Fig2showsan
exampleofabistableswitchbasedonthisidea.Itmanagestoachievetherequirednonlinearity
withonlyelementaryreactionsatthecostofarelativelymorecomplicatednetwork.Notethat
itisnotclearinrealitywhichratelawdoesapplyinmolecularregulatorynetworks.Itmightbe
neitherMichaelis-Mentennorelementaryreactionratelaws.However,multisitephosphoryla-
tionbasedonelementaryreactionratelawsdoesprovideaconvenientmechanismthatcanbe
easilymodeledeitherdeterministicallyorstochastically.Thisideawasoriginallyproposedby
Kapuyetal.[6],whoproposedthatmultisitephosphorylationsequencesmaybemodeledby
elementaryreactionratelawsandaresuitableforbothdeterministicandstochasticsimula-
tions.LaterBariketal.[7]developedamodelofyeastcellcycleregulationbasedonmultiple
site(de)phosphorylationreactionsandelementaryreactionratelaws.Theirmodelwascom-
paredwithexperimentalmeasurements[8].However,thismodelingtechniquepresentschal-
lengesforcurrentstochasticsimulationalgorithms.Mostreactionfiringsareforthe(de)
phosphorylationreactions,whichtypicallyresultinfrequentbut“back-and-forth”changes
thatcanceleachother’seffect.Theoverallstochasticcharacteristicsofthemodelismostly
drivenbyotherslowreactions.ItwouldbeidealifonecouldreducethesysteminFig2tothe
simplerversioninFig1.Studyingreductionconditionforthesystem(4)willhelpusachieve
thisgoal.
Anotherbiochemicalexampleofsystem(4)comesfromonedimensionalreaction-diffusion
system,wherediffusionisusuallymuchfasterthanreaction.Givenaparticularpartition(with
Fig2.AbistableswitchconvertedfromthemodelinFig1basedonmultiplesitephosphorylation.Cdh1isassumedtohave11phosphorylation
levels,asdescribedasCdh1(theunphosphorylatedform)andCdh1PthroughCdh1P .
10
doi:10.1371/journal.pone.0133295.g002
PLOSONE|DOI:10.1371/journal.pone.0133295 August11,2015 3/19
TheAbridgmentandRelaxationTimeforaLinearMulti-ScaleModel
nlinearbins),thediffusioncanbeformulatedasachainreaction.Supposethediffusionrates
areuniformandthebinsareequallyspaced,ifwecanreducesystem(4)tothesimplesystem
(1),thenwebasicallyreduceaspatiallyinhomogeneoussystemtoaspatiallyhomogeneoussys-
tem.Thusinthiscasetheconditionforthereductionofsystem(4)isalsorelatedtothecondi-
tionforthefamous“well-stirred”assumptiongiveninGillespie’spioneeringpaper[1,2].
Inthispaperwepresentanerroranalysisofabridgmentbasedoncomparisonofthefirst
exittime.Althoughouranalysisfocusesonthechainreactionsystem(4),theresultcanbe
extendedtomorecomplicatedlinearsystems.DetailsoftheanalysisisgivenintheMethods
sectionandwedemonstrateonespecificapplicationonbistableswitchandoscillationmodels
intheResultssection.
Methods
Background
Inthissection,webrieflyreviewthesimulationmethods:SSAandSQSSA/ssSSA.
SSA. SupposethesysteminvolvesNmolecularspecies{S ,...,S }.Thestatevectoris
1 N
denotedbyX(t)=(X (t),...,X (t)),whereX(t)isthenumberofmoleculesofspeciesS at
1 N i i
timet.Mreactionchannels{R ,...,R }areinvolvedinthesystem.Assumethatthesystemis
1 M
wellstirredandisinthermalequilibrium.ThedynamicsofreactionchannelR ischaracterized
j
bythepropensityfunctiona andthestatechangevectorν =(ν ,...,ν ):a(x)dtgivesthe
j j 1j N,j j
probabilitythatoneR reactionwilloccurinthenextinfinitesimaltimeinterval[t,t+dt),and
j
ν givesthechangeintheS moleculepopulationinducedbyoneR reaction.
ij i j
ThedynamicsofthesystemcanbesimulatedbytheSSA.WithX(t)=x,let
P
a ðxÞ¼ M aðxÞ.Ineachstep,theSSAgeneratestworandomnumbersr andr inU(0,1),
0 j¼1 j 1 2
theuniformdistributionintheinterval(0,1).Thetimeforthenextreactiontooccurisgiven
byt+τ,whereτisgivenby
(cid:2) (cid:3)
1 1
t¼ log : ð6Þ
a ðxÞ r
0 1
Theindexjforthenextreactionisgivenbythesmallestintegersatisfying
Xj
aðxÞ>r a ðxÞ: ð7Þ
l 2 0
l¼1
ThesystemstatesareupdatedbyX(t+τ)=x+ν.Thesimulationproceedstothenextoccurring
j
time,untilitreachesthefinaltimeorotherstopcriteria.
SQSSA/ssSSA. Thestochasticquasi–steady–stateapproximation(SQSSA)[9]andslow–
scaleSSA(ssSSA)[10,11]wereproposedtoreducethesizeofastochasticbiologicalmodel
andtoimprovethestochasticsimulationefficiency.ThessSSAisbasedonthepartialequilib-
rium(PE)assumption,whiletheSQSSAisbasedonthequasisteadystate(QSS)assumption.
ThePEandQSSassumptionsarebothimportantmultiscalefeaturesinbiochemicalsystems.
PEreferstothesituationwheresomereactionsfiremuchfasterthanothersandthecorre-
spondingsubsystemreachesapartialequilibriumstate[10].QSSreferstothesituationwhere
somestatevariablesfluctuateveryquicklyaroundtheirquasisteadystates[9].Ifasubsystemis
atPE,thenallitsinvolvedspeciesareinQSS.
ToimplementtheSQSSAorssSSAmethod,asystemispartitionedintofastandslowsub-
systems,andthestatevariablesarepartitionedintofastandslowvariables.Thefastsubsystem
isassumedatanequilibriumstate,whilethestatevariablesinthefastsubsystemareinQSS.
ThenthesimulationprocedureisverysimilartotheoriginalSSAprocedure,exceptthatin
PLOSONE|DOI:10.1371/journal.pone.0133295 August11,2015 4/19
TheAbridgmentandRelaxationTimeforaLinearMulti-ScaleModel
eachsimulationstep,asetofalgebraicequationsissolvedtofindtheQSSofallfastvariables.
Thenthesefastvariablesareconsidered(temporary)constantparametersintheslowsubsys-
temandtheSSAprocedureisappliedonlytheslowsubsystem.Becausethenumberofreaction
firingsintheslowsubsystemismuchlessthantheoneforthefastsubsystem,highsimulation
efficiencyisachieved.However,ifthetimescaledifferenceisnotlargeenough,theQSSA
methodmayleadtolargeerrors.Ontheotherhand,eveniftimescaleseparationholds,SQSSA
mayalsoleadtosignificanterrors.Thomasetal.studiedthisprobleminthework[12].
AnalysisfortheChainReactionSystemwithOneParticle
Erroranalysisforabridgmentinthestochasticregimeisasubtlebusiness[3].Thejustification
inSQSSA[9]andssSSA[10,11]werebothbasedondistributiondistanceofensemblescol-
lectedwithmultiplerunsoftheoriginalSSAandthecorrespondingapproximationmethods.
Butthatisnotgoodenough,inThomaset.al[13],itwasshownthat“reducedmasterequation
approachcanoverestimatethevarianceofthefluctuationsbyasmuchas*30%.Thusamore
rigorouserroranalysisisneededwhendealingwithabridgmentinthestochasticregime.In
Gillespieset.al’swork[3],anovelerroranalysismethodwasproposedbasedonthefirstexit
timeanalysis.Tostudytheconditionsothatreaction(2)isanaccuratereductionofsystem(1)
withgivenpopulationsofS andS ,thedistributionofthefirstexittime,thetimewhenthe
1 2
firstS isgenerated,wasanalyzedforbothEqs(1)and(2).Fortheabridgmenttobevalid,the
3
twodistributionsshouldbeveryclose.Inthesimplesystem(2),thefirstexittimehasanexpo-
nentialdistribution.Thustheabridgmentisaccurateonlywhenthedistributionofthefirstexit
timeforsystem(1)isclosetoanexponentialdistributionwithameanvalueclosetotheone
fromthesimplesystem(2).Otherwise,productionofS in(2)significantlydiffersfromthe
3
onein(1)andasaresultabridgmentcannotbeapplied.Conditions(3)arethusderivedbased
onthisanalysis[3].
Inthispaper,wefollowastrategysimilaraswhatwasutilizedinGillespieetal.[3]and
studythegeneralcase(4).Inordertostudythechainreactionsystem(4),westartwithasim-
plecase.SupposethereisonlyoneparticlestartingatstateS .Wewillstudythedistribution
1
functionoftheexittimeT,thefirsttimethatthisparticlebecomesS .Followingtheanalysis
n+1
ofTheorem1inGillespieet.al[3],ifthesystem(4)canbereducedto(5),thedistributionfor
Tshouldfollowanexponentialdistribution,i.e.
cdfðT;tÞ¼ProbðT <tÞ¼1(cid:4)e(cid:4)ct; ð8Þ
orequivalently
ProbðT (cid:5)tÞ¼e(cid:4)ct: ð9Þ
Denote
pðtÞ¼ProbðXðtÞ¼1Þ; for i¼1;...;nþ1; ð10Þ
i i
andletP(t)=(p (t),...,p (t))T.Notethatforthisoneparticle,ittakesoneandonlyonestate
1 n
atanygiventime.Thuswealwayshave
Xnþ1
pðtÞ¼1: ð11Þ
i
i¼1
WhenX (t)=1,itimpliesthattheexittimeT<t.Thuswehave
n+1
X
n
ProbðT <tÞ¼p ðtÞ¼1(cid:4) pðtÞ;
nþ1 i
i¼1
PLOSONE|DOI:10.1371/journal.pone.0133295 August11,2015 5/19
TheAbridgmentandRelaxationTimeforaLinearMulti-ScaleModel
orequivalently,
X
n
ProbðT (cid:5)tÞ¼ pðtÞ¼1TPðtÞ; ð12Þ
i
i¼1
where1=(1,...,1)T.Therefore,inordertoobtaintheanalyticdistributionfunctionofTwe
needtostudytheequationforP(t).WecanwritedownthechemicalmasterequationforPas:
dP
¼AP; ð13Þ
dt
where
2 3
(cid:4)f b 0 0 (cid:3)(cid:3)(cid:3) 0
1 1
6 7
6 7
6 f (cid:4)ðb þf Þ b 0 (cid:3)(cid:3)(cid:3) 0 7
6 1 1 2 2 7
6 7
6 . 7
6 . 7
6 0 f (cid:4)ðb þf Þ b . 0 7
A¼66 2 2 3 3 77; ð14Þ
6 . . . . 7
6 0 .. .. .. .. 0 7
6 7
6 7
6 7
6 0 (cid:3)(cid:3)(cid:3) 0 f (cid:4)ðb þf Þ b 7
4 n(cid:4)2 n(cid:4)2 n(cid:4)1 n(cid:4)1 5
0 (cid:3)(cid:3)(cid:3) 0 0 f (cid:4)ðb þf Þ
n(cid:4)1 n(cid:4)1 n
withtheinitialconditionP(0)=(1,0,...,0)T.
Ahasnnegativeeigenvalues[14,15].Wedenotethemasλ (cid:6)λ (cid:6)...(cid:6)λ <0,andthe
1 2 n
correspondingeigenvectorsasξ ,...,ξ .ForthelinearODE(13),wecanformulateitsanalytic
P 1 n
solutionasPðtÞ¼ ni¼1ciePlitxi,whereci’saredeterminedbytheinitialstate.Thus
ProbðT (cid:5)tÞ¼1TPðtÞ¼ ni¼1ci1Txielit,whichcanbewellapproximatedbyanexponential
functionasshownin(9)ifthefollowingconditionissatisfied:
jl j(cid:7)jlj; for all 1(cid:6)i(cid:6)n(cid:4)1: ð15Þ
n i
Moreover,whenthecondition(15)issatisfied,thereactionrateforthereducedsystem(5)isc
=jλ j.
n
Condition(15)isasufficientconditionfortheabridgmenttobevalid.Buthowisthecondi-
tion(15)relatedtothefast/slowpartitioningoftheoriginalsystem(4)?Nextwewillreveala
connectionbetweenthecondition(15)andtherelaxationtimeofthefastsubsystem.
FastSubsystemanditsRelaxationTime
Weassumethatreactionsinsystem(4)areallfastexceptthelaststep,andthenconsiderthe
fastsubsystem:
S Ðf1 S Ðf2 (cid:3)(cid:3)(cid:3)fÐn(cid:4)1S ; ð16Þ
1 2 n
b1 b2 bn(cid:4)1
withonlyoneparticleinthesystem.Forthissubsystem,wedenoteitsstatevariablesas
X^ ðtÞ;i¼1;...;n.Let
i
p^ðtÞ¼ProbðX^ ðtÞ¼1Þ; for i¼1;...;n; ð17Þ
i i
PLOSONE|DOI:10.1371/journal.pone.0133295 August11,2015 6/19
TheAbridgmentandRelaxationTimeforaLinearMulti-ScaleModel
andP^ðtÞ¼ðp^ ðtÞ;...;p^ ðtÞÞT.WehaveacorrespondingCMEforP^:
1 n
dP^ ¼BP^; ð18Þ
dt
where
2 3
(cid:4)f b 0 0 (cid:3)(cid:3)(cid:3) 0
1 1
6 7
6 7
6 f (cid:4)ðb þf Þ b 0 (cid:3)(cid:3)(cid:3) 0 7
6 1 1 2 2 7
6 7
6 . 7
6 . 7
6 0 f (cid:4)ðb þf Þ b . 0 7
B¼66 2 2 3 3 77: ð19Þ
6 . . . . 7
6 0 .. .. .. .. 0 7
6 7
6 7
6 7
64 0 (cid:3)(cid:3)(cid:3) 0 fn(cid:4)2 (cid:4)ðbn(cid:4)2þfn(cid:4)1Þ bn(cid:4)1 75
0 (cid:3)(cid:3)(cid:3) 0 0 f (cid:4)ðb Þ
n(cid:4)1 n(cid:4)1
SimilartomatrixA,Bhasn−1negativeeigenvaluesandonezeroeigenvalue[14,15].There-
fore,wecandenoteB’seigenvaluesasl^ (cid:6)l^ (cid:6)...(cid:6)l^ <l^ ¼0,andthecorresponding
1 2 n(cid:4)1 n
eigenvectorsasη ,...,η .Define
1 n
1 1
Trelax ¼1(cid:6)mi(cid:6)anx(cid:4)1jl^j¼jl^ j: ð20Þ
i n(cid:4)1
T iscalledtherelaxationtimeofthefastsubsystem(16).Itisanimportantcharacteristic
relax
forthesubsystem(16).Therelaxationtimerepresentsthetimescaleforasystemto“relax”to
itssteadystate.Italsocanbeviewedasthetimescaleafterwhichthefastsubsystemwill“for-
get”apreviousperturbation.Toseethis,westudythesolutionoftheEq(18),givenby
Xn Xn(cid:4)1
P^ ¼ ^cel^itZ ¼^c Z þ ^cel^itZ: ð21Þ
i i n n i i
1 1
Whenjl^t j(cid:2)1for1(cid:6)i(cid:6)n−1,P^ !P^(cid:8) ¼^c Z .Inotherwords,whent(cid:2)T ,P^reaches
i n n relax
itsequilibriumstateP^(cid:8).
Toseetheconnectionbetweentherelaxationtimeandtheabridgment,letE=A−B.E’sele-
mentsareallzeroexceptthelastdiagonalelement−f .Ecanbeconsideredasarankoneper-
n
turbationmatrix.Thusfortheeigenvaluesλ’sandl^’s,wehave(page101inWilkinson[16])
i i
0(cid:6)l^ (cid:4)l (cid:6)f ; for i¼1;...;n: ð22Þ
i i n
Thusjλ j(cid:6)f andjl j(cid:5)jl^ jfori=1,...,n−1.Socondition(15)willbevalidifthefollow-
n n i i
ingconditionissatisfied:
fn (cid:7)1(cid:6)mi(cid:6)inn(cid:4)1jl^ij¼jl^n(cid:4)1j: ð23Þ
Thecondition(23)hasaspecialfeature.Ontheleftside,f isthereactionratefortheslowreac-
n
^
tion,whileontherightside,jl jisacharacteristicforthefastsubsystem.Wecanalso
n(cid:4)1
PLOSONE|DOI:10.1371/journal.pone.0133295 August11,2015 7/19
TheAbridgmentandRelaxationTimeforaLinearMulti-ScaleModel
rewritecondition(23)asitsequivalentcondition:
1 1
Trelax ¼jl^ j(cid:7)f ; ð24Þ
n(cid:4)1 n
where 1 isthemeanfiringtimefortheslowreaction.
fn
Basedonthepreviousanalysis,wecanconcludethatifthecondition(23)issatisfied,the
abridgmentfromsystem(4)tosystem(5)isvalidwiththereactionratec=jλ j.
n
ReactionRatefortheReducedSystem
Theabovesectionrevealstheconnectionbetweentheabridgmentandtherelaxationtime.But
fromapracticalpointofview,westillneedtocalculatethereactionratecforthereducedsys-
tem(5).Moreover,weneedtocalculateorestimateT .
relax
First,weneedtosolvetheequilibriumstateP^(cid:8)ofthesystem(16),whichsatisfies
BP^(cid:8) ¼0: ð25Þ
Bissingular.(25)doesnotuniquelygiveasolutionfortheequilibriumstate.Weneedtorealize
thatforsystem(16),
X
n
p^ ¼1: ð26Þ
i
i¼1
^
Weapply(26)andsubstitutep toobtainanonsingularmatrix.From(26)wehave
n
Xn(cid:4)1
p^ ¼1(cid:4) p^: ð27Þ
n i
i¼1
LetP^denotethevectorfp^ ;...;p^ g.Wethenhavetheequation
1 n(cid:4)1
dP^ ¼B~P^þy; ð28Þ
dt
where
2 3
(cid:4)f b 0 0 (cid:3)(cid:3)(cid:3) 0
1 1
6 7
6 f (cid:4)ðb þf Þ b 0 (cid:3)(cid:3)(cid:3) 0 7
6 1 1 2 2 7
6 7
6 .. 7
6 0 f (cid:4)ðb þf Þ b . 0 7
~ 6 2 2 3 3 7
B ¼6 7
6 .. .. .. .. .. .. 7
6 . . . . . . 7
6 7
6 7
4 0 (cid:3)(cid:3)(cid:3) 0 f (cid:4)ðb þf Þ b 5
n(cid:4)3 n(cid:4)3 n(cid:4)2 n(cid:4)2
(cid:4)b (cid:3)(cid:3)(cid:3) (cid:4)b (cid:4)b f (cid:4)b (cid:4)ðb þf þb Þ
n(cid:4)1 n(cid:4)1 n(cid:4)1 n(cid:4)2 n(cid:4)1 n(cid:4)2 n(cid:4)1 n(cid:4)1
andy=(0,...,0,bn−1)T.OnecanverifythatB~isnonsingularandthesolutionβoftheequilib-
riumequation
B~P^þy¼0 ð29Þ
givestheuniquestableequilibriumstateforthesystem(16).
Asshowninsection,whenthecondition(15)issatisfied,thereactionrateforthereduced
system(5)isc=λ .Buthowdowecalculateλ ?Letβdenotethesolutionfor(29)andγ=
n n
PLOSONE|DOI:10.1371/journal.pone.0133295 August11,2015 8/19
TheAbridgmentandRelaxationTimeforaLinearMulti-ScaleModel
−f (1−1Tβ).Whencondition(23)issatisfied,λ canbewellapproximatedbyγintermsof
n n
BauerandFiketheorem[17]andCorollary2.2inEisenstatandIpsen[18](seeS1Appendix)
witharelativeerror
jl (cid:4)gj
n ¼Oðf T Þ:
jl j n relax
n
GeneralCase
Nowwecanloosentheconditionthatthereisonlyoneparticleinthesystem(4).Suppose
therearetotallymparticles.TheneachparticlehasitsownexittimeT,j=1,...,m.Theearli-
j
estexittimeT istheminimumofallexittimes.T isofconcernbecauseitrepresentthefirst
0 0
timeoneoftheparticleswillexitfromthesystemandmayleadtosignificanteventthatisnot
includedinthesystem.Toformulateit,wehave
T0 ¼1m(cid:6)j(cid:6)inmTj: ð30Þ
Wehave
Y
m
ProbðT0 (cid:5)tÞ¼ProbðTj (cid:5)t; for all 1(cid:6)j(cid:6)mÞ¼ ProbðTj (cid:5)tÞ: ð31Þ
j¼1
Whencondition(23)issatisfied,asshowninprevioussection,foreachj,thecorresponding
exittimeT wellapproximatesanexponentialdistributionwiththesamerateλ .Thuswehave
j n
ProbðT (cid:5)tÞ(cid:9)elnt; ð32Þ
j
and
Y
m
ProbðT0 (cid:5)tÞ¼ ProbðTj (cid:5)tÞ(cid:9)emlnt; ð33Þ
j¼1
Thusforsystem(4)withmparticles,whencondition(23)issatisfied,itcanbereducedtosys-
tem(5)withreactionratec=mjλ j.
n
Remark.
(cid:129) Adoesnothavetobeaconstantmatrix.Theaboveanalysisisbasedonalinearchainreac-
tionsystemoffastvariablesbutnotnecessaryofslowvariables.Betweenslowreactionfir-
ings,slowvariablesdonotchange.Thereactionratesmaybenonlinearfunctionsofslow
variables.Inthatcase,theeigenvaluesofthematrixAwillchangewithslowreactionfirings.
Thiswillbethecaseinournumericalexperiments.
(cid:129) Theanalysiscanbeeasilyextendedtoasystemthatisslightlydifferentfrom(4).Forexam-
ple,onecanchangethesystemto
S Ðf1 S Ðf2 (cid:3)(cid:3)(cid:3)fÐn(cid:4)1S ;
1 2 n
b1 b2 bn(cid:4)1
and
S!fn S ;
i nþ1
whereicouldbeanyintegerbetween1andn.Itiseasytoverifythatouranalysisisstillvalid.
Theonlydifferenceiswheretheextradiagonaltermf islocated.Furthermore,wewillshow
n
PLOSONE|DOI:10.1371/journal.pone.0133295 August11,2015 9/19
TheAbridgmentandRelaxationTimeforaLinearMulti-ScaleModel
inS2Appendixthatouranalysiscanalsobeextendedtosystemslike
S1ÐE S2ÐE (cid:3)(cid:3)(cid:3)ÐE Sn;
k0 k0 k0
ð34Þ
;(cid:4)!E k!uSu;;
wherethechainreactionsarefastandthedegradationreactionforEisslow.HereE,asan
enzyme,doesnotchangeitsvalueinthechainreaction,andeachindividualforwardreaction
hastheform
S þE(cid:4)!S þE; for i¼1;(cid:3)(cid:3)(cid:3);n(cid:4)1;
i iþ1
withelementaryreactionratelaws,whileS’sdonotchangevaluesintheenzymicreactions
i
EþS(cid:4)!S;
i i
either.
Results
Inthissectionwetesttheaccuracyofabridgmentwithabistableswitchmodelconsistingofa
chainreactionsystem.Numericalexperimentsontwooscillationmodelsbasedonthisbistable
switchmodelarealsopresented.Thesesimulationswereperformedona3.0GHzIntelCore
Linuxworkstation.
BistableSwitchModel
ThisisthebistableswitchmodelinFig2.Werewriteitas
Cdh1P Ðk(cid:8)Clb2 Cdh1P Ðk(cid:8)Clb2 (cid:3)(cid:3)(cid:3) Ðk(cid:8)Clb2 Cdh1P ;
0 1 9
k(cid:8)kp k(cid:8)kp k(cid:8)kp
;!ks Clb2; ð35Þ
P
9
Clb2(cid:4)ka i!¼0ð10(cid:4)iÞCdh1Pi ;:
wherek andk areconstantscorrespondingtoClb2’ssynthesisanddegradation.kisascalar
s a
thatcontrolsthereactionscaleofthechainsystem.Theseenzymicreactionsallfollowelemen-
taryreactionratelaws.Notethatouranalysiscanbeappliedtothissystem(seeS2Appendix)
althoughtheformatappearsdifferentthan(4).Wewilltesttheaccuracywhenthissystemis
reducedto
; !ksClb2;
ð36Þ
Clb2!Cdh1(cid:8)l ;;
whereCdh1isthetotalpopulationofCdh1P andλiscalculatedaccordingly.Byintuition,
i
whenkislarge,allreactionsinthelinearchainsystemarefast,thequasisteadystateassump-
tionisapplicableandasaresultthereductionisvalid.Bycontrast,ifkissosmallthatpropen-
sitiesofchainreactionsarecomparabletothetwoslowreactions,greaterrorswillbe
introduced.
PLOSONE|DOI:10.1371/journal.pone.0133295 August11,2015 10/19