A variational approach to the stochastic aspects of cellular signal transduction

YuehengLan

DepartmentofChemistry,UniversityofNorthCarolinaatChapelHill,NC27599-3290

PeterG.Wolynes

DepartmentofChemistry&Biochemistry,UniversityofCaliforniaatSanDiego,9500GilmanDr.,LaJolla,CA92093-0371

GareginA.Papoian∗

arXiv:q-bio/0607025v1 [q-bio.QM] 18 Jul 2006DepartmentofChemistry,UniversityofNorthCarolinaatChapelHill,NC27599-3290

(Dated:February5,2008)

Cellularsignalingnetworkshaveevolvedtocopewithintrinsicfluctuations,comingfromthesmallnumbersofconstituents,andtheenvironmentalnoise.Stochasticchemicalkineticsequationsgovernthewaybiochemicalnetworksprocessnoisysignals.Theessentialdifficultyassociatedwiththemasterequationapproachtosolvingthestochasticchemicalkineticsproblemistheenormousnumberofordinarydifferentialequationsinvolved.Inthiswork,weshowhowtoachievetremendousreductioninthedimensionalityofspecificreactioncascadedynamicsbysolvingvariationallyanequivalentquantumfieldtheoreticformulationofstochasticchemicalkinetics.Thepresentformulationavoidscumbersomecommutatorcomputationsinthederivationofevolutionequations,makingmoretransparentthephysicalsignificanceofthevariationalmethod.Weproposenoveltime-dependentbasisfunctionswhichworkwelloverawiderangeofrateparameters.Weapplythenewbasisfunctionstodescribestochasticsignalinginseveralenzymaticcascadesandcomparetheresultssoobtainedwiththosefromalternativesolutiontechniques.Thevariationalansatzgivesprobabilitydistributionsthatagreewellwiththeexactones,evenwhenfluctuationsarelargeanddiscretenessandnonlinearityareimportant.AnumericalimplementationofourtechniqueismanyordersofmagnitudemoreefficientcomputationallycomparedwiththetraditionalMonteCarlosimulationalgorithmsortheLangevinsimulations.

Keywords:StochasticProcesses,NonlinearChemicalKinetics,VariationalApproach,QuantumFieldTheory,SignalTrans-duction,DiscreteNoise,StrongFluctuations,MasterEquation

I.INTRODUCTION

Thelifeofthecellisregulatedbyintricatechainsofchemicalreactions1.Thewholecellmaybeviewedasacomputingdevicewhereinformationisreceived,relayedandprocessed2.Signaltransductioncascadesbasedonproteininteractionsreg-ulatecellmovement,metabolismanddivision1,3.Sincecellsaremesoscopicobjects,understandingtheroleoftheintrinsicfluctuationsofthebiochemicalreactionsaswellasenvironmentalfluctuationsisafundamentalpartofunderstandingsignalingdynamics4,5,6,7,8,9,10,11,12,13,14,15.Inthisregard,thewell-organizedbehaviorofcells,whichemergesasaresultofbiochemicalreactiondynamicsinvolvinghundredsofcross-linkedsignalingpathways,isremarkable16,17,18,19,20,21,22.Theproblemofhowsignalscanbepreciselydetected,smoothlytransducedandreliablyprocessedundernoisyconditionsisaresearchtopicofgreatcurrentinterest,that,inturn,shouldleadtodeeperunderstandingoftheoriginsofthecell’sfunctionalresponses23,24.Further-more,thesestudiescanhelptounravelthedesignprinciplesforvarioussignalingpathways,leading,eventually,tobetterwaystocontrolandefficientlyinterferewithcellularactivity,aswouldbeneededtocorrectthebehaviorofdiseasedcells18,25.

Theroleofnoiseingeneregulatorynetworkshasbeenidentifiedasakeyissueandhasbeenintensivelystudiedinrecentyears10,11,12,26,27,28,29,30,31.Linearizationofthenoisemaybeacceptableifthedynamicsnearsteadystatesisbeingstudied10,26,31.Whenproteinnumbersarelargeand,thus,thecontinuousapproximationisvalid,time-dependentdistributionshavebeende-terminedusingtheLangevinorFokker-Planckequations6,32,33.Toaccountforthediscretenessinthelinearizedequations,thegeneratingfunctionapproachhasalsobeenused10,26.Avariationaltreatmentofsteadystatestabilityandswitchinginnonlinear,discretegeneregulatoryprocesseshasbeenreported29,30.

Incytosolicsignaltransductionprocesses,incontrasttogenetranscriptionwhichinvolvesauniqueDNAmolecule,allthereactingspeciesarepresentinmultiplecopiesandparticipateinunary,binaryorperhapsevenhigherorderreactions.Noisecouldbemultiplicative34,35andthelineardescriptioneasilybreaksdown.Moreover,cellularreactionsusuallytakeplacehetero-geneouslyinspaceThelocalizationandcompartmentalizationofproteinorganellesrequirediffusiveoractivetransportationofreactingmoleculesfromoneregiontoanother.Spatialcoordinationcombinedwithtemporalcoordinationgeneratescoherent,yetcomplexspatiotemporalpatterns18,19,20,21,22,36,37,38.

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

dt

(m,n)=µ[−mnP(m,n)+m(n+1)P(m,n+1)]

+λ[−(N−n)p(m,n)+(N−n+1)P(m,n−1)]

+g[−P(m,n)+P(m−1,n)]+k[−mP(m,n)+(m+1)P(m+1,n)],

(1)

whereNisthetotalnumberofAandA∗.InEq.(1),thefirsttwotermsdescribetheA−A∗reactionandtheresttheR−R∗reaction.Thissimple2-stepcascadeiscommonlyfoundembeddedintheonsetofareactionpathwayofmanyimportantsignalingcascades55,56.Ifalargenumberofinactivereceptors,R,arepresenttherateofconversiondependsonthearrivaltimesoftheexternalcueandthereactionbecomesPoissonian.Weassumethatthisisthecaseinallthefollowingcalculations.IftheR→R∗reactionistheusualbirth-deathproblem,ourformalismstillapplieswithonlyminorchanges.Themasterequation(1)actuallycontainsinfinitelymanycoupledODEs.

B.TheQFTformulation

Thedifferential-differenceequations,suchasEq.(1),arewellrepresentedintheQFTformulationbyintroducingcreationandannihilationoperatorsa,a†andstates|n󰀎29,30,47,48,49,50,51.Inanalogytoquantummechanics,theoperatorssatisfythecom-mutationrelationthat

[a,a†]=1.

Asusual,the“vacuumstate”|0󰀎anditsconjugate󰀐0|aredefinedtosatisfy

󰀐0|a†=a|0󰀎=0,

󰀐0|0󰀎=1.

Otherstatesarebuiltupfromthevacuumstate,suchasthen-particlestate|n󰀎

|n󰀎=a†n|0󰀎.

3

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

=Ω|Ψ󰀎.(2)dt

TheoriginallargesetsofODEsarenowcompactifiedintojustoneequation.Appliedtothe2-stepcascade(Fig.1),Eq.(2)ischaracterizedbyfollowingoperatorΩ,

whereb†,barethecreationandannihilationoperatorsassociatedwithR∗anda†,awithA.Inthiscase,

󰀁|Ψ󰀎=P(m,n)|m,n󰀎,

m,n

Ω=(1−a†)(µb†ba−λN+λa†a)+g(b†−1)+k(b−b†b),

(3)

where

Eq.(3)isreadilyverifiedbysubstitutingintoEq.(2)andcomparingthecoefficientsofeach(m,n)-particlestate.Incontrastto

ordinaryquantummechanics,theoperatorΩisnon-Hermitian,sotheinnerproductsbetweenthestatesarenotconserved.Thiswasthereasontointroduceearliertheharvestingstate.Nevertheless,manyQFTtechniquesmaybeprofitablyapplied,albeitwithsomemodifications29,30,48,50,57,58,59.Wewillnotdiscussthoseandinsteadwilltranslatetheabovefieldtheoreticformulationtothefamiliardifferentialequationlanguage.

C.Differentialoperatorsandthegeneratingfunction

|m,n󰀎=a†mb†n|0󰀎.

Inthefieldtheoreticform,thecomputationsarecarriedoutbycommutatormanipulationsthatsometimesareawkward.Fortunately,itturnsoutthatwemayconverttheoperatorequation(2)intoapartialdifferentialequation(PDE).Toaccomplishthat,weexploretheanalogybetweena,a†andd/dx,x.Notonlydotheyhavethesamecommutator

[a,a†]=1⇐⇒

[d

a|0󰀎=0⇐⇒

dx

dxd

x=1;

xn=nxn−1;

ama†n|0󰀎=

n!

dx

)mxn=

n!

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

amf(a†)|0󰀎⇐⇒(

d

dx

f(x);

∂t

=(1−x)(µy

∂2

∂x

)Ψ+g(y−1)Ψ−k(y−1)

∂Ψ

|x=1,∂x∂2Ψ2

󰀐n󰀎=(

∂x

)|x=1,(5)

where|x=1meansevaluationatx=1.Therefore,themomentsareobtainedwhenthegeneratingfunctionisexpandedatx=1,whiletheprobabilitydistributionisobtainedfromtheTaylorcoefficientswhenthegeneratingfunctionΨisexpandedatx=0.

D.Thevariationalsolution

IntheQFTformulationofthestochasticprocesses,avariationalprinciplemaybederivedwhichisequivalenttotheevolutionequation(2).ThisprincipleindicatesthatthephysicalsolutionofEq.(2)isgivenbythestationarypointsofthefollowingfunctional29,30,53

󰀂∞

H[ΦL,ΦR]=dt󰀐ΦL|∂t−Ω|ΦR󰀎,(6)

0

whereΦLandΦRarearbitraryquantumstatesunderquitegeneralconstraintsconsistentwiththepositivityofprobabilitiesand

thefixedboundaryconditions.Inpractice,wetakeafinite-dimensionalsubsetoftheinfinite-dimensionalfunctionspaceandapplythevariationalprincipleinthissubspacetogetclosedequationsthatmaybesubsequentlysolvedbysimplenumericalcalculation.Iftheessentialqualitativepropertiesofthesystemareknown,goodapproximationsoftheoriginalproblemcanbeachievedthroughaninformedchoiceoftime-dependentbasisfunctionsthatdefinetherelevantsubsetinthefunctionspace.BecauseΩisnotHermitian,therightandlefteigenvectorsaredifferent.Tocharacterizethesystem,we,therefore,needtwosetsofvectorsΦLandΦR.ThestationaryvariationconditionforΦLrestorestheoriginalequation(2)andthatforΦRdefinesanequationsatisfiedbyΦL.Ifweviewtheoperator∂t−Ωasalargematrixparameterizedbyt,theΦLandΦRgeneratedbythestationaryvariationconditioncorrespondtoitssingularvectors61,62andtheextremumvaluesofEq.(6)arethesingularvalues.Physically,fromtheSchr¨odingerpicturepointofview,ΦRistheevolvingquantumstateandΦLrepresentsthemeasurablequantitiesinwhichweareinterested.Eq.(6)servestofindthemostsignificantstateandphysicalobservables.EyinkoriginallyappliedthisvariationalprincipletoFokker-Planckequations54.Subsequently,SasaiandWolynesusedthisvariationalapproachinthefieldtheoreticformandobtainedmomentsinatoggle-switchgeneregulatoryproblem29.Inthispaper,weshowhowthevariationalprinciplemaybeapplied,instead,tothegeneratingfunctions.Weintroducenovelbasisfunctionstoobtainthetime-dependentprobabilitydistributionsinsignaltransductioncascades.Anothernoveltyofthepresentformulationisouravoidance

5

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

dxi

ThefunctionalHsimplybecomes

H=

󰀂

∞0

.(9)

dtΦL(∂t−Ω)ΦR|x=1,

(10)

Inthenewpicture,wehavemuchsimplermathematicaloperations,e.g.,thevariationalprinciplebecomessimplyafunctionextremizationcondition

δH

dxi

(∂t−Ω)ΦR|x=1=0,,fori=1,2,...,m.

(12)

Theevolutionofthegeneratingfunctionshouldalwaysconservethetotalprobability.AsinEq.(4),thetotalprobabilityΨ(1,1)doesnotchangewithtime.TheproperchoiceofΦRshouldalsoguaranteethisinvariance,satisfying(∂t−Ω)ΦR|x=1=0.Now,Eq.(12)tellsusthatthehigherderivativesoftheexpression(∂t−Ω)ΦRevaluatedatx=1arealsozero.Therefore,inthelimitofm→∞,Eq.(12)leadstothePDE∂tΨR=ΩΨR.Forfinitem,thisPDEisapproximatelysatisfiedintheneighborhoodofx=1.

III.

NUMERICALAPPLICATIONS

Inthispart,wediscusstheimplementationofthePDEversionofthevariationalmethodandapplyittoseveralsimple,yetimportantenzymaticcascades.Beforeproceedingtotheindividualexamples,weemphasizeourmotivationforselectingthetime-dependentbasisfunctions.Wealsobrieflydiscussseveralalternativemethodsalsousedtosolvethemasterequation.

A.

Computationaldetails

ItisreasonabletorequirethefollowingconstraintsontherightbasisfunctionΦR(x,y):(1)thetotalprobabilityshouldbeequalto1,i.e.,ΦR(1,1)=1;

(2)theprobabilityshouldbepositive,i.e.,thecoefficientsoftheTaylorexpansionofΦRaround(x,y)=(0,0)shouldbenonnegative;

(3)thetimerateoftheunknownfunctions,f˙i(t),shouldbeobtainablebysolvingEq.(12)derivedfromthevariationalprinciple.Inthefollowing,weintroducetwosetsofbasisfunctions.Onesetissimplebutisoflimitedapplicability,whiletheotherisinamorecomplexintegralformandcanbeappliedverygenerally.

6

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

λ+µm(t)

followingansatz,

k(1

+

µm(t)

−e−kt)istheaveragenumberofR∗attimet.Wemakeuseofthisspecificfunctionalformandtrythe

ΦR=(1+f2(t)(y−1))Ne(x−1)f1(t),

(18)

whichresultsinthefollowing2DODEs

f˙1=g−kf1

f˙2=λ(1−f2)−µf1f2.

(19)

Theseequationshaveaparticularlysimplephysicalexplanation-theycorrespondtothedeterministicchemicalkineticsequa-tionssincef1andN∗f2areequaltotheaveragenumbersofR∗andA,respectively.Butnow,wemayobtainprobability

2

distributionsthroughEq.(18).Forexample,thevarianceofAcanbeeasilycalculatedasσ2=󰀐n2󰀎−󰀐n󰀎2=f2(t)−f2(t).

∗

TheseODEscanbesolvedexactlyandweshowinFig.2theprobabilitydistributionofAatt=30fortwosetsofparametervalues.Alsoshowninthefigureareresultsobtainedfromcalculationsusingmoretraditionaltechniques.Thefirstsetofreactionrateparameterswerechosenasg=10,k=5,µ=0.02,λ=0.15,withtheinitialconditions(NR,NR∗,NA,NA∗)=(20,0,5,0).SincetheR−R∗reactionismuchfasterthantheA−A∗reaction,oneexpectsEq.(17)tobeagoodapproximation40.Indeed,inFig.2(a),thevariationalansatzEq.(18)leadstoaresultthatoverlapssignificantlybetterwiththeexactGillespiecalculation,comparedwiththeresultsfromΩ-expansionandLangevinequation.TheΩ-expansionresultturnsouttobemoreconcentratedthantheexactresult,whiletheLangevinequationdoesnotworkwellneartheleftboundary,shiftingtheaveragetotheright.

Forotherparametervalues,aslongastheR−R∗reactionisfast,theansatzEq.(18)worksfineasexpected40.However,ifthefirstreactionisconsiderablyslowerthanthesecondone,thisansatzbecomeslessuseful,asshowninFig.2(b)for

7

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

√

√

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

∂2

∂t

=(1−z)(−µ2y+(1−y)(µx

∂2

√

∂z

)Ψ

)Ψ+g(x−1)Ψ−k(x−1)

∂Ψ

∂y

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

∂z∂w

−λ3N3+(λ3w+µ3N2)

∂

√

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction

Y.Lan,P.G.Wolynes,&G.A.PapoianAvariationalapproachtostochasticsignaltransduction