Injektivitás és multistacionaritás, összetett kémiai reakciókban

Az SR-gráf

Példaként egy nemtriviális (tehát a program algoritmusai le kell fussanak rajta) reakciólista (hogy ne kelljen mindent mindig újra beírni):
R={{“a+m”,”<->”,”x”},{“b+n”,”<->”,”y”},{“2*a+n”,”<-“,”y”},{“b+x”,”<->”,”z”},{“r+m”,”<-“,”z”}}

R={{"a+m","<->","x"},{"b+n","<->","y"},{"2*a+n","<-","y"},{"b+x","<->","z"},{"r+m","<-","z"}}

A következő függvények a reakciólistából elkészítik a gráf éleit. Éllistás formátumban tároltam, mert először nem tudtam pontosan milyen formában érdemes kezelnem. Amíg a feladatot nem oldom meg hiánytalanul, megtartom ezt a formátumot, hátha még szükségem lesz rá. Megj.: Minden függvény után van egy parancs ami kipróbálja azt, így könnyen lehet látni mit is csinál pontosan.

edges[r_]:=Module[{R1=StringSplit[#,"+"]&/@r,R2={},s,e},{Do[AppendTo[R2,StringSplit[#,"*"]&/@R1[[i]]],{i,1,Length[R1]}],s=species[R2],e=Partition[s,1],For[l=0,l<Length[R2],l++;e=addedge[R2[[l]],l,e,r]]};{s,StringJoin/@r,e}](*az eredménylista első eleme: az anyagcsúcsok, második eleme: a reakciócsúcsok, harmadik eleme: az élek listája.*)

species[r_]:= (*Kigyűjti a szereplő anyagokat a reakciókból.*)Module[{s={},rr=Flatten[r,2]},For[j=0,j<Length[rr],j++;If[rr[[j]][[1]]=="->"||rr[[j]][[1]]=="<-"||rr[[j]][[1]]=="<->",,If[Length[rr[[j]]]>1,s=Union[s,{rr[[j]][[2]]}],s=Union[s,{rr[[j]][[1]]}]]]];s]
species[{{{{"d"}},{{"->"}},{{"e"},{"3","f"}}},{{{"2","a"},{"b"}},{{"<->"}},{{"d"}}}}]

addedge[r_,l_,s_,re_]:=Module[{rr=Flatten[r,1],a,pos,irany=1,coef=1,e=s}, (*Az r reakciót hozzáadja az s listához (a megfelelő anyaghoz fűzi élként), ahol az r-et a re reakciólista l-edik elemeként kaptuk.*)For[k=0,k<Length[rr],k++;If[rr[[k]][[1]]=="<-"||rr[[k]][[1]]=="->"||rr[[k]][[1]]=="<->",irany=3,If[Length [rr[[k]]]>1,{coef=ToExpression[rr[[k]][[1]]],a=rr[[k]][[2]],pos=Position[e,a,2][[1]][[1]],AppendTo[e[[pos]],{StringJoin[re[[l]]],re[[l]][[irany]],coef}]},{coef=1,a=rr[[k]][[1]],pos=Position[e,a,2][[1]][[1]],AppendTo[e[[pos]],{StringJoin[re[[l]]],re[[l]][[irany]],coef}]}]]];e]

edgelist=edges[R]

Graphlist={};For[m=0,m<Length[edgelist[[3]]],m++;For[mm=1,mm<Length[edgelist[[3]][[m]]],mm++;AppendTo[Graphlist,edgelist[[3]][[m]][[1]]\[UndirectedEdge]edgelist[[3]][[m]][[mm]][[1]]]]]

Itt létrehozok egy gráf objektumot, ettől a ponttól kezdve tudom használni a beépített parancsokat. A DFS-re a körkereséshez lesz szükség. Itt sajnos elakadtam, mivel egy gráfban az összes kört megkeresni nehéz feladat, ezért mindenképpen fel kell használnom a gráf speciális tulajdonságait, illetve, hogy sok kört meg se kell vizsgálni. Még nem sikerült rájönnöm hogyan tudnám ezt hatékonyan elvégezni, de remélem a félév végéig ez sikerülni fog.

G=Graph[Graphlist]

Egy rövid kis segédfüggvény, majd két elővizsgálat: Ha megfigyeljük a tételt, az (i) feltétel teljesül ha minden él sztoichiometrikus együtthatója 1, illetve (ii) teljesül ha minden anyagcsúcs foka legfeljebb 3 (ezek le vannak írva a cikkben is).
Ezért az alábbi két kis program ellenőrzi ezen feltételek teljesülését.

not1[n_]:=True
not1[n_]:=False /; NumberQ[n] && n!=1

Prechk1[ed_]:=Module[{check=True},check=Not[MemberQ[not1/@Flatten[ed],False]];check];

Prechk2[g_,l_]:=Not[MemberQ[#>3&/@(VertexDegree[g,#]&/@l),True]]

Prechk1[edgelist]

Prechk2[G,edgelist[[1]]]

Sajnos itt elakadtam, és nem sikerült rájönnöm hogyan lehetne hatékonyan megoldani a feltételek ellenőrzését. Bővebben kifejtem a mellékelt latex dokumentumban.

A polinomfüggvény injektivitásának vizsgálata:

Az alábbiakban rendszeresen fogok utalni olyan eredményekre, amiket részletesen leírtam a mellékelt latex dokumentumban.

A továbbiakban egy másik módszert programozok le. A bemenetről feltesszük, hogy a reakciók minden esetben mint n dimenziós vektorok vannak megadva, ahol n jelöli az anyagfajták számát. Ezt az n-et könnyen ki lehet olvasni a reakciólistából, de az egyszerűség kedvéért inkább úgy írtam meg a programot hogy feltesszük, hogy előre ismert.
Pl: S={A,B,C,D,X,Y,Z,W};
R={{{1,1,0,0,0,0,0,0},{0,0,1,0,0,0,0,0}},{{0,0,1,0,0,0,0,0},{1,1,0,0,0,0,0,0}}, {{0,0,0,0,1,0,0,0},{2,0,0,1,0,0,0,0}},{{2,0,0,1,0,0,0,0},{0,0,0,0,1,0,0,0}}, {{2,0,0,1,0,0,0,0},{0,0,0,0,0,1,0,0}},{{0,0,0,0,0,1,0,0},{2,0,0,1,0,0,0,0}}, {{0,0,0,1,0,0,0,0},{0,0,1,0,0,0,0,1}},{{0,0,1,0,0,0,0,1},{0,0,0,1,0,0,0,0}}, {{0,1,0,1,0,0,0,0},{0,0,0,0,0,0,1,0}},{{0,0,0,0,0,0,1,0},{0,1,0,1,0,0,0,0}}, {{1,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0}},{{0,1,0,0,0,0,0,0},{0,0,0,0,0,0,0,0}}, {{0,0,1,0,0,0,0,0},{0,0,0,0,0,0,0,0}},{{0,0,0,1,0,0,0,0},{0,0,0,0,0,0,0,0}}, {{0,0,0,0,1,0,0,0},{0,0,0,0,0,0,0,0}},{{0,0,0,0,0,1,0,0},{0,0,0,0,0,0,0,0}}, {{0,0,0,0,0,0,1,0},{0,0,0,0,0,0,0,0}},{{0,0,0,0,0,0,0,1},{0,0,0,0,0,0,0,0}}}

Ez az első program True-val tér vissza ha a reakcióhálózat injektív, ha pedig nem, akkor lefut a *firstsubcheck alparancs, ami ha False-t ad, akár csak egyszer is (ezért catch-throw-t használtam), az azt jelenti, hogy a reakcióhálózat biztosan multistacionárius. Ha nem ad False-t, akkor le kell futtatni a *secondchecket*, és az még szolgálhat eredménnyel. Ezt a 2. Megjegyzésben részleteztem.

firstcheck[r_,n_]:=Module[{lr=Length[r],result=True,results,subs1,subs2},{subs1=Map[List[#,1]&,Subsets[Range[lr],{n}],{2}],subs2=Map[List[#,2]&,Subsets[Range[lr],{n}],{2}],If[Catch[results=Parallelize[Map[If[Det[Extract[r,subs1[[#]]]]*Det[Extract[r,subs1[[#]]]-Extract[r,subs2[[#]]]]<0,If[Not[firstsubcheck[Transpose[Extract[r,subs1[[#]]]-Extract[r,subs2[[#]]]],n]],Throw[False],1],0]&,Range[Binomial[lr,n]]]]]==False,result=False],If[Count[results,1]>0 && result≠False,result="Indecisive"]};result]

firstsubcheck[y_,n_]:=Module[{result=True,t},{tab=Table[t,{n}],If[FindInstance[Inner[GreaterEqual,y.tab,Table[1,{n}],And]&&Inner[GreaterEqual,tab,Table[1,{n}],And],tab]≠{},result=False]};result]

Timing[firstcheck[{{{1,1,0,0,0,0,0,0},{0,0,1,0,0,0,0,0}},{{0,0,1,0,0,0,0,0},{1,1,0,0,0,0,0,0}},{{0,0,0,0,1,0,0,0},{2,0,0,1,0,0,0,0}},{{2,0,0,1,0,0,0,0},{0,0,0,0,1,0,0,0}},{{2,0,0,1,0,0,0,0},{0,0,0,0,0,1,0,0}},{{0,0,0,0,0,1,0,0},{2,0,0,1,0,0,0,0}},{{0,0,0,1,0,0,0,0},{0,0,1,0,0,0,0,1}},{{0,0,1,0,0,0,0,1},{0,0,0,1,0,0,0,0}},{{0,1,0,1,0,0,0,0},{0,0,0,0,0,0,1,0}},{{0,0,0,0,0,0,1,0},{0,1,0,1,0,0,0,0}},{{1,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0}},{{0,1,0,0,0,0,0,0},{0,0,0,0,0,0,0,0}},{{0,0,1,0,0,0,0,0},{0,0,0,0,0,0,0,0}},{{0,0,0,1,0,0,0,0},{0,0,0,0,0,0,0,0}},{{0,0,0,0,1,0,0,0},{0,0,0,0,0,0,0,0}},{{0,0,0,0,0,1,0,0},{0,0,0,0,0,0,0,0}},{{0,0,0,0,0,0,1,0},{0,0,0,0,0,0,0,0}},{{0,0,0,0,0,0,0,1},{0,0,0,0,0,0,0,0}}},8]]

Megjegyzés: összehasonlítás képpen, parallelize nélkül az algoritmus több mint tízszer ennyi ideig futott. bemásolom ide, tesztelhetőség kedvéért:

regifirstcheck[r_,n_]:=Module[{lr=Length[r],result=True,subs1,subs2,i=1},{subs1=Map[List[#,1]&,Subsets[Range[lr],{n}],{2}],subs2=Map[List[#,2]&,Subsets[Range[lr],{n}],{2}],While[result ≠ False && i≤Binomial[lr,n],If[Det[Extract[r,subs1[[i]]]]*Det[Extract[r,subs1[[i]]]-Extract[r,subs2[[i]]]]<0,result=regifirstsubcheck[Transpose[Extract[r,subs1[[i]]]-Extract[r,subs2[[i]]]],n]];i++]};result]
regifirstsubcheck[y_,n_]:=Module[{result="Indecisive",t},{tab=Table[t,{n}],If[FindInstance[Inner[GreaterEqual,y.tab,Table[1,{n}],And]&&Inner[GreaterEqual,tab,Table[1,{n}],And],tab]≠{},result=False]};result]

A fenti példára tehát választ ad az első algoritmus. A következőkben olyan példát tekintünk, amire a fenti algoritmus nem ad választ, de a második, erősebb igen.
S2= {A,B}; R2={{{2,1},{3,0}},{{3,0},{2,1}},{{1,0},{0,0}},{{0,1},{0,0}}}

firstcheck[{{{2,1},{3,0}},{{3,0},{2,1}},{{1,0},{0,0}},{{0,1},{0,0}}},2]

secondcheck[r_,n_]:=Module[{result="Indecisive",μ,t,lr=Length[r]},{μ=Array[t,lr],T[mu_,δ_]:=Sum[mu[[i]] (r[[i,1]].δ) (r[[i,1]]-r[[i,2]]),{i,lr}],f[mu_]:=Det[ToMatrix[T,mu,n]],If[FindInstance[f[μ]<0 && Inner[Greater,μ,Table[0,{lr}],And] && Inner[Greater,Sum[μ[[i]](r[[i,1]]-r[[i,2]]),{i,lr}],Table[0,{n}],And],μ]≠{},result=False]};result]

ToMatrix[f_,mu_,n_]:=Table[f[mu,UnitVector[n,i]],{i,n}]

secondcheck[{{{2,1},{3,0}},{{3,0},{2,1}},{{1,0},{0,0}},{{0,1},{0,0}}},2]

És végül a kettőt egyesítem egy programban, ami a fentiekhez hasonlóan értelmezhető eredményt ad, tehát True ha injektív (és akkor egyértelmű is), False ha biztosan multistacionárius, és Indecisive, ha nem injektív ugyan, de nem is biztosan multistacionárius.

MultEqQ[r_,n_]:=Module[{result=firstcheck[r,n]},{If[result=="Indecisive",result=secondcheck[r,n]]};result]

MultEqQ[{{{2,1},{3,0}},{{3,0},{2,1}},{{1,0},{0,0}},{{0,1},{0,0}}},2]

Egy további példa amivel a program leellenőrízhető (ezek a C-F cikkből vannak, az 1.1-es táblázatban találhatók):

A táblázatban (iv): S4={A,B,C,P,Q} ; R4={{{1,1,0,0,0},{0,0,0,1,0}},{{0,0,0,1,0},{1,1,0,0,0}},{{0,1,1,0,0},{0,0,0,0,1}}, {{0,0,0,0,1},{0,1,1,0,0}},{{0,0,1,0,0},{1,0,0,0,0}},{{1,0,0,0,0},{0,0,1,0,0}},{{1,0,0,0,0},{0,0,0,0,0}}, {{0,1,0,0,0},{0,0,0,0,0}},{{0,0,1,0,0},{0,0,0,0,0}},{{0,0,0,1,0},{0,0,0,0,0}}, {{0,0,0,0,1},{0,0,0,0,0}}}

MultEqQ[{{{1,1,0,0,0},{0,0,0,1,0}},{{0,0,0,1,0},{1,1,0,0,0}},{{0,1,1,0,0},{0,0,0,0,1}},{{0,0,0,0,1},{0,1,1,0,0}},{{0,0,1,0,0},{1,0,0,0,0}},{{1,0,0,0,0},{0,0,1,0,0}},{{1,0,0,0,0},{0,0,0,0,0}},{{0,1,0,0,0},{0,0,0,0,0}},{{0,0,1,0,0},{0,0,0,0,0}},{{0,0,0,1,0},{0,0,0,0,0}},{{0,0,0,0,1},{0,0,0,0,0}}},5]

Gheorghe Craciun, Martin Feinberg,Multiple equilibria in complex chemical reaction networks: I. The injectivity property, SIAM J. Appl. Math.,Vol. 65, No. 5, pp. 1526-1546,2005; Link

Gheorghe Craciun, Martin Feinberg,Multiple equilibria in complex chemical reaction networks: II. The species-reaction graph, SIAM J. Appl. Math.,Vol. 66, No. 4, pp. 1321-1338,2006; Link