University of Groningen Enabling Darwinian evolution in chemical replicators Mattia, Elio

(1)

University of Groningen

Enabling Darwinian evolution in chemical replicators

Mattia, Elio

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from

it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date:

2017

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Mattia, E. (2017). Enabling Darwinian evolution in chemical replicators. Rijksuniversiteit Groningen.

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the

author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately

and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the

number of authors shown on this cover page is limited to 10 maximum.

(2)

149

149 Ena

bl

ing

D

ar

w

ini

an e

vo

lu

tio

n i

n c

he

m

ic

al

re

pl

ic

at

or

s

A.

APPENDIX:

SOURCE CODE

(3)

150 Appe

nd

ix

A. S

ou

rc

e c

od

e

_{In this appendix, source code is provided for the various models developed for this}

research and presented in this thesis.

The sections are organized as follows:

A.1.

C++ deterministic model, exponential replication (Chapter 2)

A.2. Matlab deterministic model, far-from-equilibrium replication (Chapter

3, Section 3.2)

A.2.1. Main code (single and multiple parameter sets)

A.2.2. Kinetic equation set

A.2.3. Flux equation set

A.2.4. 3D plot toolset

A.2.5. Parameters panels

A.3. C++ stochastic model (Chapter 3, Section 3.3)

Section A.2 presents the toolset used in the Matlab model in detail, as indicated

above.

A.1. C++ DETERMINISTIC MODEL,

EXPONENTIAL REPLICATION (CHAPTER 2)

// :: KinSim :: //

// Author: Elio Mattia // Date: February 2011 (v1) // Required files: none //

/*

The second version has a different model with respect to KinSim: the first version bypasses solution hexamers and allows direct nucleation of a stack of two hexamers from solution.

The 2.0 version doesn't allow for nucleation from trimers and tetramers and considers hexamers in solution, which form from trimers and tetramers. Only solution hexamers can nucleate, from hexamers, then two growth mechanisms are considered: sequestration of hexamers from solution (which is exponential at first, then can become linear, rate-limited by the rate of solution-hexamers formation from the library) and catalyzed synthesis of hexamers, just from trimers as a first model, which is exponential if fibre breakage is considered.

VERSION HISTORY

2.1 - Flag "steady" added: it allows, if on, to maintain a steady state of monomer to its initial concentration, instead of it being consumed

2.2 - Added backward reactions from 6 towards 3 and 4

3.0 - Allows calculation of P thanks to outputting log(d[R]/dt) against log[R], the slope is P, the order of reaction. "steady" must be on.

3.1 - Introduced "steady6" flag, keeps steady the amount of **solution** hexamers. Keeping steady monomer doesn't allow calculation of P

3.2 - Introduced "steady3" flag, useful for catalyzed elongation simulations to calculate P

3.3 - Introduced fibre length distribution histograms calculation and printout 3.3a - Not a full version, set maxlen to 100 and thr1 to 10 instead of maxlen/2, longer simulations with different breakage probability than the standard one

3.4 - Introduced 'nucleationon' flag, to (de)activate nucleation from solution hexamers; IMPORTANT: if deactivated, automatically seeds the initial system with a constant distribution of 0.1 to all fibre lengths

3.5 - Allowing, thanks to "fibrefeedback", for fibre concentration to influence breakage rate. The law considers a term as independent from fibre concentration and another one as dependent from it.

*/ //--- #include <math.h> #include <stdio.h> #include <stdlib.h> #include <string.h> #include <time.h>

(4)

151 Appe

nd

ix

A. S

ou

rce

co

de

//=============================== :: Typedef Struct :: ================================== typedefstruct configuration {

float s1; float s2; float s3; float s4; float s6; float*stacks; } configuration; //============================== :: Global Constants :: ================================= constint v1=3; constint v2=5; constchar vl=' '; constint maxlen=100; constint maxlength=maxlen+1; constlong maxtime=10000; constfloat monomerinitconc=1; constfloat trimerinitconc=2e2; constfloat hexamerinitconc=1e1; constint steady=0;

constint steady3=1; constint steady6=1; constfloat rate12=1e-1; constfloat rate23=1e-1; constfloat rate24=3e-1; constfloat rate36=3e-1; constfloat rate46=3e-1; constfloat rate63=1e-5; constfloat rate64=1e-5; constfloat nuclprob=1e-6; constfloat sequrate=0; constfloat elonrate=3e-1; constint nucleationon=1; constint breakageon=0; constint fibrefeedback=0; constint fft=10; constint thr1=maxlen/2; constint thr2=maxlen;

//================================== :: Prototypes :: ===================================

void simulate(float*,configuration *); void cfgcopy(configuration *,configuration *); void generatebp(float*);

configuration *create(); void destroy(configuration *); void createpath(char*);

void write(long,FILE*,configuration *,float*,float*,float*); void parse(FILE*,char*);

//======================================================================================= //======================================================================================= //===================================== :: MAIN :: ====================================== //======================================================================================= //======================================================================================= int main(){

//==================== :: structures and structured lists :: ========================

configuration *cfg; //============================ :: local variables :: ================================ long t; char path[250]; char txtpath[254]; FILE*output; float bp[maxlength]; int i; float r61,r62,r63; //================================= :: main code :: ================================== generatebp(bp); createpath(path); strcpy(txtpath,path); strcat(txtpath,".txt"); output=fopen(txtpath,"w");

fprintf(output,"dynamics_%d%c%d N%d T%d C%.0e 12_%.0e 23_%.0e 24_%.0e 36_%.0e 46_%.0e 63_%.0e 64_%.0e NP%.0e SR%.0e ER%.0e break%d

steady%d\n",v1,'.',v2,maxlen,maxtime,monomerinitconc,rate12,rate23,rate24,rate36,rate46,rate63,rate6 4,nuclprob,sequrate,elonrate,breakageon,steady);

fprintf(output,"dynamics_%d%c%d\tN%d\tT%d\tC%.0e\t12_%.0e\t23_%.0e\t24_%.0e\t36_%.0e\t46_%.0e\t63_%. 0e\t64_%.0e\tNP%.0e\tSR%.0e\tER%.0e\tbreak%d\tsteady%d\n\n",v1,'.',v2,maxlen,maxtime,monomerinitconc ,rate12,rate23,rate24,rate36,rate46,rate63,rate64,nuclprob,sequrate,elonrate,breakageon,steady);

(5)

152 Appe

nd

ix

A. S

ou

rc

e c

od

e

fprintf(output,"%c\t%d\t%d\t%d\t%d\t%d%c%c%c%c\t%d%c%c%c%c\t\t%d%c%c%c%c\t",'t',1,2,3,4,6,' ','s','o','l',6,' ','f','i','b',6,' ','t','o','t');

for(i=1;i<maxlength;i++) fprintf(output,"\t%d",i); fprintf(output,"\t\tlog[R]\tlog(d[R]/dt)\t"); for(i=1;i<maxlength;i++) fprintf(output,"\t%d",i); fprintf(output,"\n\n"); r61=r62=r63=0; cfg=create(); cfg->s1=monomerinitconc; cfg->s3=trimerinitconc; cfg->s6=hexamerinitconc; t=0;

write(t,output,cfg,&r61,&r62,&r63);

for(++t;t<=maxtime;t++){ simulate(bp,cfg);

//Move information relevant to calculation of p

r61=r62; r62=r63;

//Print results

write(t,output,cfg,&r61,&r62,&r63); }

destroy(cfg);

fclose(output); output=fopen(txtpath,"r"); parse(output,path); fclose(output); return0; } //======================================================================================= //==================================== :: SIMULATE :: =================================== //=======================================================================================

void simulate(float*bp,configuration *cfg){ configuration *next; int i; float totstacks; float breakprob; next=create();

//UPDATES: rates proportional to reactants concentrations and kinetic constants

//update element (increase), then update also the source (decrease),

//consider stoichiometry: 6mM of 1 equals 3mM of 2, 2mM of 3, 1,5 mM of 4, 1 mM of 6

//also, 1 mM of 6_2 equals 1mM of fibre ends, 1 mM of 6_3 equals 0,66 mM of fibre ends, 1 mM of 6_4 equals 0,5 mM of fibre ends, etc. (so, it's 2/n)

//Update 2 (from 1) if(cfg->s1>=rate12*(cfg->s1)*(cfg->s1)){ next->s2=cfg->s2+rate12*(cfg->s1)*(cfg->s1)/2; if(steady){ next->s1=cfg->s1; }else{ next->s1=cfg->s1-rate12*(cfg->s1)*(cfg->s1); } }else{ next->s2=cfg->s2+cfg->s1/2; if(steady){ next->s1=cfg->s1; }else{ next->s1=0; } }

//Update 3 (from 2 and 1 combined)

if((next->s2>=rate23*(cfg->s2)*(cfg->s1))&&(next->s1>=rate23*(cfg->s2)*(cfg->s1))){ if(steady3){

next->s3=cfg->s3; }else{

next->s3=cfg->s3+rate23*(cfg->s2)*(cfg->s1); }

next->s2=next->s2-rate23*(cfg->s2)*(cfg->s1); if(steady){

next->s1=next->s1; }else{

next->s1=next->s1-rate23*(cfg->s2)*(cfg->s1); }

}else{

(6)

153 Appe

nd

ix

A. S

ou

rce

co

de

if(steady3){ next->s3=cfg->s3; }else{ next->s3=cfg->s3+next->s2; } if(steady){ next->s1=next->s1; }else{

next->s1=next->s1-next->s2; } next->s2=0; }else{ if(steady3){ next->s3=cfg->s3; }else{ next->s3=cfg->s3+next->s1; }

next->s2=next->s2-next->s1; if(steady){ next->s1=next->s1; }else{ next->s1=0; } } } //Update 4 (from 2) if(next->s2>=rate24*(cfg->s2)*(cfg->s2)){ next->s4=cfg->s4+rate24*(cfg->s2)*(cfg->s2)/2; next->s2=next->s2-rate24*(cfg->s2)*(cfg->s2); }else{

next->s4=cfg->s4+next->s2/2; next->s2=0;

}

//Update 6 (from 3, and from 4 and 2 combined)

if(next->s3>=rate36*(cfg->s3)*(cfg->s3)){ if(steady6){ next->s6=cfg->s6; }else{ next->s6=cfg->s6+rate36*(cfg->s3)*(cfg->s3)/2; } if(steady3){ next->s3=next->s3; }else{

next->s3=next->s3-rate36*(cfg->s3)*(cfg->s3); } }else{ if(steady6){ next->s6=cfg->s6; }else{ next->s6=cfg->s6+next->s3/2; } if(steady3){ next->s3=next->s3; }else{ next->s3=0; } }

if((next->s4>=rate46*(cfg->s4)*(cfg->s2))&&(next->s2>=rate46*(cfg->s4)*(cfg->s2))){ if(steady6){

next->s6=next->s6+rate46*(cfg->s4)*(cfg->s2); }

next->s4=next->s4-rate46*(cfg->s4)*(cfg->s2); next->s2=next->s2-rate46*(cfg->s4)*(cfg->s2); }else{

if(next->s2>=next->s4){ if(steady6){ next->s6=next->s6; }else{

next->s6=next->s6+next->s4; }

next->s2=next->s2-next->s4; next->s4=0;

}else{ if(steady6){ next->s6=next->s6; }else{

next->s6=next->s6+next->s2; }

next->s4=next->s4-next->s2; next->s2=0; } } //Update 3 (from 6) if(next->s6>=rate63*(cfg->s6)){ if(steady3){

(7)

154 Appe

nd

ix

A. S

ou

rc

e c

od

e

next->s3=next->s3+rate63*(cfg->s6)*2; }

if(steady6){ next->s6=next->s6; }else{

next->s6=next->s6-rate63*(cfg->s6); }

}else{ if(steady3){ next->s3=next->s3; }else{

next->s3=next->s3+next->s6*2; } if(steady6){ next->s6=next->s6; }else{ next->s6=0; } } //Update 4 (from 6) if(next->s6>=rate64*(cfg->s6)){ next->s4=next->s4+rate64*(cfg->s6)*3/2; if(steady6){

next->s6=next->s6-rate64*(cfg->s6); }

}else{

next->s4=next->s4+next->s6*3/2; if(steady6){ next->s6=next->s6; }else{ next->s6=0; } }

//Update conc of all the lengths towards +1 growth: for all other lengths except 2,

next->stacks[maxlength-1]=cfg->stacks[maxlength-1];//working always on next avoids update problems

for(i=maxlength-1; i>=3; i--){

//consider ===> SEQUESTRATION <=== with sequestration rate constant, FIBRE END CONC ([conc of (n-1)], times 2, divided by n-1), conc of 6

next->stacks[i-1]=cfg->stacks[i-1];

if((next->stacks[i-1]>=sequrate*cfg->stacks[i-1]*2/(i-1)*(cfg->s6)*(i-1))&&(next ->s3>=sequrate*cfg->stacks[i-1]*2/(i-1)*(cfg->s6))){

next->stacks[i]=next->stacks[i]+sequrate*cfg->stacks[i-1]*2/(i-1)*(cfg->s6)*i; //hexamer moles become hexamers-on-stack, and take with them (i-1) times moles of hexamers involved in fibres (i-1) hexamers long,

//to form fibres i hexamers long

next->stacks[i-1]=next->stacks[i-1]-sequrate*cfg->stacks[i-1]*2/(i-1)*(cfg->s6)*(i-1); if(steady6){

next->s6=next->s6-sequrate*cfg->stacks[i-1]*2/(i-1)*(cfg->s6); }

}else{

if(next->s6>=next->stacks[i-1]/(i-1)){

next->stacks[i]=next->stacks[i]+next->stacks[i-1]*i/(i-1); if(steady6){

next->s6=next->s6-next->stacks[i-1]/(i-1); }

next->stacks[i-1]=0; }else{

next->stacks[i]=next->stacks[i]+next->s6*i; next->stacks[i-1]=next->stacks[i-1]-next->s6*(i-1); if(steady6){ next->s6=next->s6; }else{ next->s6=0; } } }

//and ===> CATALYZED ELONGATION <=== with elongation rate constant, fibre end conc (same as above), conc of 3

if((next->stacks[i-1]>=elonrate*cfg->stacks[i-1]*2/(i-1)*(cfg->s3)*(cfg->s3)/2*(i -1))&&(next->s3>=elonrate*cfg->stacks[i-1]*2/(i-1)*(cfg->s3)*(cfg->s3))){

next->stacks[i]=next->stacks[i]+elonrate*cfg->stacks[i-1]*2/(i-1)*(cfg->s3)*(cfg ->s3)/2*i;

//half of trimer moles become hexamers, and take with them (i-1) times moles of hexamers involved in fibres (i-1) hexamers long,

//to form fibres i hexamers long

next->stacks[i-1]=next->stacks[i-1]-elonrate*cfg->stacks[i-1]*2/(i-1)*(cfg->s3)*(cfg ->s3)/2*(i-1);

(8)

155 Appe

nd

ix

A. S

ou

rce

co

de

if(steady3){ next->s3=next->s3; }else{

next->s3=next->s3-elonrate*cfg->stacks[i-1]*2/(i-1)*(cfg->s3)*(cfg->s3); }

}else{

if(next->s3>=next->stacks[i-1]*2/(i-1)){

next->stacks[i]=next->stacks[i]+next->stacks[i-1]*i/(i-1); if(steady3){

next->s3=next->s3-next->stacks[i-1]*2/(i-1); }

next->stacks[i-1]=0; }else{

next->stacks[i]=next->stacks[i]+next->s3/2*i; next->stacks[i-1]=next->stacks[i-1]-next->s3/2*(i-1); if(steady3){ next->s3=next->s3; }else{ next->s3=0; } } } }

//for length 2, consider nucleation rate and conc of 6

if(nucleationon){

if(next->s6>=nuclprob*(cfg->s6)*(cfg->s6)){

next->stacks[2]=next->stacks[2]+nuclprob*(cfg->s6)*(cfg->s6); if(steady6){

next->s6=next->s6-nuclprob*(cfg->s6)*(cfg->s6); }

}else{

next->stacks[2]=next->stacks[2]+next->s6; if(steady6){ next->s6=next->s6; }else{ next->s6=0; } } }

//calculate total fibre concentration

totstacks=0;

for(i=1; i<maxlength; i++){ totstacks+=cfg->stacks[i]; }

//Update conc of lower lengths upon breakage of higher lengths: rate proportional to probability (high number approx: we

//don't consider probability distributions, but just rates), update higher length, then lower length(s)

if(breakageon){

for(i=2; i<maxlength; i++){ if(fibrefeedback){

breakprob=(bp[i]*(1+totstacks/fft)>1)?1:(bp[i]*(1+totstacks/100000)); }else{

breakprob=bp[i]; }

if(!(i%2)){

next->stacks[i/2]=next->stacks[i/2]+breakprob*next->stacks[i]; next->stacks[i]=next->stacks[i]-breakprob*next->stacks[i]; }else{

next->stacks[(int)(i/2)]=next->stacks[(int)(i/2)]+breakprob*next->stacks[i]*(i -1)/(2*i);

next->stacks[(int)(i/2)+1]=next->stacks[(int)(i/2)+1]+breakprob*next ->stacks[i]*(i+1)/(2*i);

next->stacks[i]=next->stacks[i]-breakprob*next->stacks[i]; } } } cfgcopy(cfg,next); destroy(next); } //======================================================================================= //=================================== :: CFGCOPY :: ===================================== //=======================================================================================

void cfgcopy(configuration *c1,configuration *c2){ int i;

(9)

156 Appe

nd

ix

A. S

ou

rc

e c

od

e

c1->s1=c2->s1; c1->s2=c2->s2; c1->s3=c2->s3; c1->s4=c2->s4; c1->s6=c2->s6; for(i=0; i<maxlength; i++){ c1->stacks[i]=c2->stacks[i]; } } //======================================================================================= //================================== :: GENERATEBP :: =================================== //=======================================================================================

void generatebp(float*bp){ int i;

//Linear ramp 0-1 from thr1 to thr2, zeroing before

for(i=1; i<=thr1; i++){ bp[i]=0; } for(; i<=thr2; i++){ bp[i]=(float)(i-thr1)/(float)(thr2-thr1); } } //======================================================================================= //==================================== :: CREATE :: ===================================== //======================================================================================= configuration *create(){ configuration *cfg; int i;

cfg=(configuration *)(malloc(sizeof(configuration))); cfg->s1=0;

cfg->s2=0; cfg->s3=0; cfg->s4=0; cfg->s6=0;

cfg->stacks=(float*)(malloc(maxlength*sizeof(float))); for(i=0;i<maxlength;i++){ cfg->stacks[i]=nucleationon?0:0.1; } return cfg; } //======================================================================================= //=================================== :: DESTROY :: ===================================== //=======================================================================================

void destroy(configuration *cfg) { free(cfg->stacks); free(cfg); } //======================================================================================= //================================== :: CREATEPATH :: =================================== //=======================================================================================

void createpath(char *txtpath) { FILE*lastpath;

lastpath=fopen("lastpath.txt","w");

fprintf(lastpath,"dyn%d%c%d%c N%d T%d 1C%.0e 3C%.0e 6C%.0e 12_%.0e 23_%.0e 24_%.0e 36_%.0e 46_%.0e 63_%.0e 64_%.0e NP%.0e SR%.0e ER%.0e n%d b%d ff%d fft%d 1t%d 2t%d 1s%d 3s%d

6s%d",v1,'.',v2,vl,maxlen,maxtime,monomerinitconc,trimerinitconc,hexamerinitconc,rate12,rate23,rate2 4,rate36,rate46,rate63,rate64,nuclprob,sequrate,elonrate,nucleationon,breakageon,fibrefeedback,fft,t hr1,thr2,steady,steady3,steady6);

fclose(lastpath);

lastpath=fopen("lastpath.txt","r"); fgets(txtpath,250,lastpath); fclose(lastpath); }

//======================================================================================= //==================================== :: WRITE :: ====================================== //=======================================================================================

void write(long t,FILE *output,configuration *cfg,float *r61,float *r62,float *r63) {

(10)

157 Appe

nd

ix

A. S

ou

rce

co

de

float totstacks; float tot; totstacks=0; for(i=1;i<maxlength;i++){ totstacks+=cfg->stacks[i]; } tot=cfg->s1+2*cfg->s2+3*cfg->s3+4*cfg->s4+6*cfg->s6+6*totstacks; *r63=6*totstacks; fprintf(output,"%d\t%f\t%f\t%f\t%f\t%f\t%f\t\t%f\t",t,cfg->s1,2*cfg->s2,3*cfg->s3,4*cfg ->s4,6*cfg->s6,6*totstacks,tot); for(i=1;i<maxlength;i++){

fprintf(output,"\t%f",cfg->stacks[i]); }

if(t>=2){

if((*r63!=*r61)&&(*r62!=0)){

fprintf(output,"\t\t%f\t%f\t",(float)(log((double)(*r62))),(float)(log((double)(*r63 )-(double)(*r61))));

for(i=1;i<maxlength;i++){ if(totstacks!=0){

fprintf(output,"\t%f",100/totstacks*(cfg->stacks[i])); } } } } fprintf(output,"\n"); } //======================================================================================= //===================================== :: PARSE :: ===================================== //=======================================================================================

void parse(FILE*output,char*path){ FILE*final;

char xlspath[154]; char character;

strcpy(xlspath,path); strcat(xlspath,".xls"); final=fopen(xlspath,"w");

while(!(feof(output))){ fscanf(output,"%c",&character); if(character=='.'){ fprintf(final,","); }else{

fprintf(final,"%c",character); } } fclose(final); } //======================================================================================= //======================================================================================= //=================================== :: END OF CODE :: ================================= //======================================================================================= //=======================================================================================

(11)

158 Appe

nd

ix

A. S

ou

rc

e c

od

e

A.2. MATLAB DETERMINISTIC MODEL,

FAR-FROM-EQUILIBRIUM REPLICATION

(CHAPTER 3, SECTION 3.2)

A.2.1 Main code (single and multiple parameter sets)

Single parameter set

% ODEsim 1.0 % Fibres model % Elio Mattia % 25/3/2014 %init clear close all clc %time tic

%core model constants

maxlen=100; %EVEN INTEGER for proper kB calculation!

n_ode=7+maxlen-1; n_flux=21;

%ODE options

options = odeset('AbsTol', 1e-12); options = odeset(options, 'RelTol', 1e-12); options = odeset(options, 'NonNegative', 1:n_ode);

%Plot options

LineWidth = 5; FontSize = 22;

%===PARAMETERS PANEL===

%flow start and total simulation time (1e4 minutes = 1 week)

fstart= 2*1e4; %start the flow

ttotal= 6*1e4; %stop the simulation

%concentration

vctot=3.8e-3; % LAB = 3.8 mM (variable)

vcox=0; %batch concentration of oxidizing agent at start of flow

%kinetic constants (per second/per minute)

kO2 = 1e-3; %GUESS

kred = 1.1e7; % LAB = 1.9e5/1.1e7, pH 8

kox = 69; % LAB = 1.15/69, pH 8

kquench = 900; % LAB = 15/900, pH 8

kflow = 5e-6; % LAB (variable) = 8.0e-6 M min-1 = 3.8 mM in 8 h

kflowox = 1e0*kflow; % LAB (variable)

kflowred = 1e0*kflow; % LAB (variable)

kX = 4.77e5; % LAB = 7951/4.77e5, pH 8

kXF = 0; % INACTIVE: exchange with fibres

kN = 1e-7; %GUESS (replication LAB = 0.593/35.6)

kS = 0; % INACTIVE: sequestration

kCE = 4e3; %GUESS (replication LAB = 0.593/35.6)

kBm = 1e2; %GUESS (replication LAB = 0.593/35.6)

%redox selectivities

sel_red = 0.3; % 1=only destroy hexamers

sel_ox = 0.2; % 1=only produce hexamers %===/PARAMETERS PANEL=== %breakage profile kB = zeros(1,maxlen); for i=(maxlen/2):maxlen kB(i)=(i-maxlen/2)/(maxlen-maxlen/2)*kBm; end %fibre length fiblen=zeros(1,maxlen); for i=1:maxlen fiblen(i)=i; end

(12)

159 Appe

nd

ix

A. S

ou

rce

co

de

%init output structures

finaly=zeros(length(fstart),length(vcox),length(kflow),n_ode); flux=zeros(n_flux,length(fstart),length(vcox),length(kflow)); flux01=zeros(length(fstart),length(vcox),length(kflow)); flux02=zeros(length(fstart),length(vcox),length(kflow)); flux03=zeros(length(fstart),length(vcox),length(kflow)); flux04=zeros(length(fstart),length(vcox),length(kflow)); flux05=zeros(length(fstart),length(vcox),length(kflow)); tElapsed=zeros(length(fstart),length(vcox),length(kflow)); %MAIN CYCLE for jctot=1:length(vctot)

disp(['ctot ',num2str(jctot),'/',num2str(length(vctot))]); disp(['---']);

disp([' ']);

for jfstart=1:length(fstart)

disp(['fstart ',num2str(jfstart),'/',num2str(length(fstart))]); for jcox=1:length(vcox)

disp([' cox ',num2str(jcox),'/',num2str(length(vcox))]); for jkflow=1:length(kflow)

tStart=tic; %time computation is for a ctot:kX:kXF:kD quadruplet (kX variation not stored)

%init time (mesh per minute)

t=linspace(0,ttotal,ttotal);

%init total concentration

ctot=vctot(jctot);

%init starting conditions

y0=zeros(n_ode,1); y0(1)=ctot; y0(7)=ctot*vcox(jcox); %init parameters par=zeros(1,12+maxlen); par(1)=maxlen; par(2)=kO2; par(3)=kred; par(4)=kox; par(5)=kquench; par(6)=kflowox(jkflow); par(7)=kflowred(jkflow); par(8)=kX; par(9)=kXF; par(10)=kN; par(11)=kS; par(12)=kCE; for k=2:maxlen par(12+k-1)=kB(k); end par(12+maxlen)=fstart(jfstart); par(12+maxlen+1)=sel_red; par(12+maxlen+2)=sel_ox; %solve ODEs [t,y]=ode15s(@gFibres5_4,t,y0,options,par); finaly(jfstart,jcox,jkflow,:)=y(length(t),:);

%compute final fluxes

flux(:,jfstart,jcox,jkflow)=fluxFibres5_4(y(length(t),:),par);

%process final fluxes

%flux01(jfstart,jcox,jkflow)=(flux(jkflow,5)+flux(jkflow,6)+flux(jkflow,8))/flux(jkflow,3);

%flux02(jfstart,jcox,jkflow)=(flux(jkflow,5)+flux(jkflow,6))/flux(jkflow,8);

%signal end of computations

tElapsed(jctot,jcox,jkflow)=toc(tStart);

disp([' kflow ',num2str(jkflow),'/',num2str(length(kflow)),': ',num2str(y(length(t),7)/ctot)]);

disp([' red34vs6 ',num2str(flux(15,jfstart,jcox,jkflow),7),' : ',num2str(flux(16,jfstart,jcox,jkflow),7)]);

disp([' ox34vs6 ',num2str(flux(17,jfstart,jcox,jkflow),7),' : ',num2str(flux(18,jfstart,jcox,jkflow),7)]);

disp([' redvsredox ',num2str(flux(19,jfstart,jcox,jkflow),7)]); disp([' replred ',num2str(flux(20,jfstart,jcox,jkflow),7)]); disp([' quench ',num2str(flux(21,jfstart,jcox,jkflow),7)]); disp(' '); disp('---'); disp([strrep(num2str(sel_red,7),'.',','),char(9),strrep(num2str(sel_ox,7),'.',',')]); disp(' '); disp([strrep(num2str(flux(20,jfstart,jcox,jkflow),7),'.',',' ),char(9),strrep(num2str(1-flux(20,jfstart,jcox,jkflow),7),'.',',')]);

(13)

160 Appe

nd

ix

A. S

ou

rc

e c

od

e

disp([strrep(num2str(flux(17,jfstart,jcox,jkflow),7),'.',','),char(9),strrep(num2str(flux(18,jfsta rt,jcox,jkflow),7),'.',',')]); disp([strrep(num2str(flux(15,jfstart,jcox,jkflow),7),'.',','),char(9),strrep(num2str(flux(16,jfsta rt,jcox,jkflow),7),'.',',')]); disp(' '); disp([strrep(num2str(flux(21,jfstart,jcox,jkflow),7),'.',',')]); disp([strrep(num2str(kflow,7),'.',',')]); disp('---'); %PLOTS

%kinetics plots (fibres and redox agents on panels 1 and 2, respectively)

figure(1)

hue=jkflow/length(kflow);

color=[(hue) 3*(hue)*(1-hue) (1-hue)]; %B->G->R (with increasing kflow)

hFig = figure(1);

set(hFig, 'Position', [50 100 1800 800]) subplot(1,2,1)

if (fstart(jfstart)>=0) && (fstart(jfstart)<=ttotal)

plot([fstart(jfstart) fstart(jfstart)],[0 1],':','Color',[0 0

0],'LineWidth',LineWidth)

end

hold on for k=1:7 if (k==1)

plot(t,y(:,k)/ctot,'Color',[1 0 0],'LineWidth',LineWidth) else

if (k==2)

if (k==5)

if (k==6) subplot(1,2,2)

plot(t,y(:,k)/ctot,'Color',[1 0 0],'LineWidth',LineWidth)

%reductant red else if (k==7) for i=1:fstart(jfstart) y(i,k)=0; end plot(t,y(:,k)/ctot,'Color',[0 0 0.5],'LineWidth',LineWidth) %oxidant dark blue

subplot(1,2,1) else

plot(t,y(:,k)/ctot,'Color',[0 0 1],'LineWidth',LineWidth) end end end end end hold on end totfib=zeros(length(t),1); fib=zeros(maxlen,1); fibdis=zeros(length(y),maxlen); for k=7+1:7+maxlen-1 totfib=totfib+y(:,k); fib(k-7)=y(length(t),k); fibdis(:,k-7)=y(:,k); end

plot(t,totfib/ctot,'Color',[0 0 0],'LineWidth',LineWidth) hold on

subplot(1,2,1)

set(gca, 'FontSize', FontSize, 'LineWidth', LineWidth) xlabel('time (minutes)')

ylabel('relative concentrations') subplot(1,2,2)

if (fstart(jfstart)>=0) && (fstart(jfstart)<=ttotal) ymax=round(1.2*(y(ttotal,7))/ctot,2,'significant');

plot([fstart(jfstart) fstart(jfstart)],[0 ymax],':','Color',[1 1

yl=ylim;

plot([fstart(jfstart) fstart(jfstart)],[0 yl(2)],':','Color',[0 0

end

set(gca, 'FontSize', FontSize, 'LineWidth', LineWidth) xlabel('time (minutes)')

ylabel('relative concentrations')

%replication constant plot (on figure 2)

dtotfib=zeros(length(t),1); krepl=zeros(length(t),1); for k=2000:20000

(14)

161 Appe

nd

ix

A. S

ou

rce

co

de

dtotfib(k,1)=totfib(k,1)-totfib(k-1,1); krepl(k,1)=dtotfib(k,1)/((y(k,3)+y(k,4))*(y(k,3)+y(k,4))*totfib(k,1)); end %figure(2) %plot(t,krepl,'Color',[1 0 0],'LineWidth',LineWidth) %hFig = figure(2); %set(hFig, 'Position', [50 100 1800 800]) %hold on disp([' ']);

disp(['krepl ',num2str(krepl(15000,1))]); disp([' ']);

%displaying ratio redox/replication

disp([' ']);

disp(['redox/replication

',num2str((flux(3)+flux(4))/(flux(7)+flux(8)+flux(9)+flux(10)),5)]); disp([' ']);

%final fibre length plot (on panel 3)

%figure(1) %subplot(1,3,3) %plot(fiblen,fib/ctot,'Color',color,'LineWidth',LineWidth) %xlabel('fibre length') %ylabel('relative amount') %axis tight %set(gca,'XTick',[0 10 20 30 40 50 60 70 80 90 100]) end end end end %figure(2) %end, time and save

disp('Ending!'); tTotal=toc; cd('D:\PhD\ODEsim\2015-2016') save('kineticsFibres_5_4.mat');

Multiple parameter set

% ODEsim 1.0 % Fibres model % Elio Mattia % 25/3/2014 %init clear close all clc %time tic

maxlen=100; %EVEN INTEGER for proper kB calculation!

n_ode=7+maxlen-1; n_flux=21;

%ODE options

options = odeset('AbsTol', 1e-12); options = odeset(options, 'RelTol', 1e-12); options = odeset(options, 'NonNegative', 1:n_ode);

%Plot options

LineWidth = 2;

%flow start and total simulation time (1e4 minutes = 1 week)

fstart= 2*1e4; %start the flow

ttotal= 6*1e4; %stop the simulation %concentration

vctot=3.8e-3; % LAB = 3.8 mM (variable)

vcox=0; %batch concentration of oxidizing agent at start of flow %kinetic constants (per second/per minute)

kO2 = 1e-3; %GUESS

kred = 1.1e7; % LAB = 1.9e5/1.1e7, pH 8

kox = 69; % LAB = 1.15/69, pH 8

kquench = 900; % LAB = 15/900, pH 8

kflow = logspace(-10,-5,11); % LAB (variable) = 8.0e-6 M min-1 = 3.8 mM in 8 h

kflowox = 1e0.*kflow; % LAB (variable)

(15)

162 Appe

nd

ix

A. S

ou

rc

e c

od

e

kX = 4.77e5; % LAB = 7951/4.77e5, pH 8

kXF = 0; % INACTIVE: exchange with fibres

kN = 1e-7; %GUESS (replication LAB = 0.593/35.6)

kS = 0; % INACTIVE: sequestration

kCE = 4e3; %GUESS (replication LAB = 0.593/35.6)

kBm = 1e2; %GUESS (replication LAB = 0.593/35.6) %redox selectivities

sel_red = linspace(0,1,11); % 1=only destroy hexamers

sel_ox = linspace(0,1,11); % 1=only produce hexamers %breakage profile kB = zeros(1,maxlen); for i=(maxlen/2):maxlen kB(i)=(i-maxlen/2)/(maxlen-maxlen/2)*kBm; end %fibre length fiblen=zeros(1,maxlen); for i=1:maxlen fiblen(i)=i; end

%init output structures

finaly=zeros(length(kflow),length(sel_red),length(sel_ox),n_ode); flux=zeros(n_flux,length(kflow),length(sel_red),length(sel_ox)); sshexa=zeros(length(kflow),length(sel_red),length(sel_ox)); sshexa=zeros(length(kflow),length(sel_red),length(sel_ox)); flux01=zeros(length(kflow),length(sel_red),length(sel_ox)); flux02=zeros(length(kflow),length(sel_red),length(sel_ox)); flux03=zeros(length(kflow),length(sel_red),length(sel_ox)); flux04=zeros(length(kflow),length(sel_red),length(sel_ox)); flux05=zeros(length(kflow),length(sel_red),length(sel_ox)); flux06=zeros(length(kflow),length(sel_red),length(sel_ox)); tElapsed=zeros(length(kflow),length(sel_red),length(sel_ox)); %MAIN CYCLE for jctot=1:length(vctot)

disp(['ctot ',num2str(jctot),'/',num2str(length(vctot))]); disp(['---']);

disp([' ']);

disp([num2str(length(kflow)),char(9),char(9),num2str(length(sel_red)),char(9),char(9),num2str(leng th(sel_ox))]);

disp(['kflow',char(9),'sel_red',char(9),'sel_ox']); for jkflow=1:length(kflow)

for jsel_red=1:length(sel_red) for jsel_ox=1:length(sel_ox)

disp([num2str(jkflow),char(9),char(9),num2str(jsel_red),char(9),char(9),num2str(jsel_ox)]); tStart=tic; %time computation is for a ctot:kX:kXF:kD quadruplet (kX variation not stored)

%init time (mesh per minute)

t=linspace(0,ttotal,ttotal);

%init total concentration

ctot=vctot(jctot);

%init starting conditions

y0=zeros(n_ode,1); y0(1)=ctot; y0(7)=ctot*vcox; %init parameters par=zeros(1,12+maxlen); par(1)=maxlen; par(2)=kO2; par(3)=kred; par(4)=kox; par(5)=kquench; par(6)=kflowox(jkflow); par(7)=kflowred(jkflow); par(8)=kX; par(9)=kXF; par(10)=kN; par(11)=kS; par(12)=kCE; for k=2:maxlen par(12+k-1)=kB(k); end par(12+maxlen)=fstart; par(12+maxlen+1)=sel_red(jsel_red); par(12+maxlen+2)=sel_ox(jsel_ox); %solve ODEs [t,y]=ode15s(@gFibres5_4,t,y0,options,par);

(16)

163 Appe

nd

ix

A. S

ou

rce

co

de

finaly(jkflow,jsel_red,jsel_ox,:)=y(length(t),:);

%compute final fluxes

flux(:,jkflow,jsel_red,jsel_ox)=fluxFibres5_4(y(length(t),:),par);

%signal end of computations

tElapsed(jkflow,jsel_red,jsel_ox)=toc(tStart); %processing totfib=zeros(length(t),1); fib=zeros(maxlen,1); fibdis=zeros(length(y),maxlen); for k=7+1:7+maxlen-1 totfib=totfib+y(:,k); fib(k-7)=y(length(t),k); fibdis(:,k-7)=y(:,k); end sshexa(jkflow,jsel_red,jsel_ox)=totfib(length(t))/ctot; %replication constant plot data

dtotfib=zeros(length(t),1); krepl=zeros(length(t),1); for k=2000:20000 dtotfib(k,1)=totfib(k,1)-totfib(k-1,1); krepl(k,1)=dtotfib(k,1)/((y(k,3)+y(k,4))*(y(k,3)+y(k,4))*totfib(k,1)); end end end end end %figure(2) %end, time and save

disp('Ending!'); tTotal=toc; cd('D:\PhD\ODEsim\2015-2016') save('kineticsFibres5_4_multiRun.mat');

A.2.2 Kinetic equation set

% ODEsim 1.0 % Fibres model % Elio Mattia % 25/3/2014

function dy=gFibres5_4(t,y,par) %#ok<INUSL> %fetching maximum fibre length

maxlen = par(1);

%calculating number of ODEs

n_ode = 7+maxlen-1; %preallocating memory kB=zeros(1,maxlen); hexafib=zeros(1,maxlen); dhexafib=zeros(1,maxlen);

%fetching kinetic constants

kO2 = par(2); kred = par(3); kox = par(4); kquench = par(5); kflowox = par(6); kflowred = par(7); kX = par(8); kXF = par(9); kN = par(10); kS = par(11); kCE = par(12); for k=2:maxlen kB(k) = par(12+k-1); end fstart = par(12+maxlen); sel_red = par(12+maxlen+1); sel_ox = par(12+maxlen+2);

%do not flow until fstart

(17)

164 Appe

nd

ix

A. S

ou

rc

e c

od

e

kred=0; kox=0; kflowox=0; kflowred=0; end

%enter glove box when starting flow: no O2 oxidation

if (t>=fstart) kO2=0; end %fetching variables mono=y(1); di=y(2); tri=y(3); tetra=y(4); hexa=y(5); red=y(6); ox=y(7); for i=2:maxlen hexafib(i)=y(7+i-1); end

%limiting oxidation rate

lim_ox=1e6;

%COMPUTING DIFFERENCES

%oxidation kinetics and exchange kinetics

if (kO2*mono>lim_ox) mono_ox=lim_ox; else mono_ox=kO2*mono; end if (kO2*di>lim_ox) di_ox=lim_ox; else di_ox=kO2*di; end

dmono= -mono_ox -kox*ox*mono;

ddi= mono_ox +kox*ox*mono di_ox di_ox kox*ox*di -kox*ox*di;

dtri= di_ox +(1-sel_ox)*kox*ox*di -kX*tri*(mono+di) +kX*tetra*(mono+di);

dtetra= di_ox +(1-sel_ox)*kox*ox*di +kX*tri*(mono+di) -kX*tetra*(mono+di);

dhexa=0;

dox= 0.5*kox*ox*mono 0.5*kox*ox*di -0.5*kox*ox*di;

%inflow and quench kinetics

dred= kflowred -kquench*red*ox; dox= dox +kflowox -kquench*red*ox;

for i=2:maxlen if (i==2)

%nucleation kinetics

dtri = dtri -kN*tri*tri;

dtetra = dtetra -kN*tetra*tetra;

dhexafib(i) = dhexafib(i) +kN*tri*tri +kN*tetra*tetra +2*sel_ox*kox*ox*di; else

%sequestration kinetics

dhexa = dhexa -kS *hexa*2/(i-1)*hexafib(i-1); dhexafib(i-1) = dhexafib(i-1) -kS*(i-1)*hexa*2/(i-1)*hexafib(i-1); dhexafib(i) = dhexafib(i) +kS*i *hexa*2/(i-1)*hexafib(i-1); %catalyzed elongation kinetics

dtri = dtri -kCE* tri*tri*2/(i-1)*hexafib(i-1); dtetra = dtetra -kCE *tetra*tetra*2/(i-1)*hexafib(i-1);

dhexafib(i-1) = dhexafib(i-1) 1)*tri*tri*2/(i-1)*hexafib(i-1) -kCE*(i-1)*tetra*tetra*2/(i-1)*hexafib(i-1);

dhexafib(i) = dhexafib(i) +kCE*i *tri*tri*2/(i-1)*hexafib(i-1) +kCE*i *tetra*tetra*2/(i-1)*hexafib(i-1);

end

if (i>=4)

%fibre breakage kinetics

if mod(i,2) %odd fibres

dhexafib(i) = dhexafib(i) -kB(i)*hexafib(i); dhexafib((i-1)/2) = dhexafib((i-1)/2) +kB(i)*hexafib(i)/2; dhexafib((i-1)/2+1) = dhexafib((i-1)/2+1) +kB(i)*hexafib(i)/2; else

%even fibres

dhexafib(i) = dhexafib(i) -kB(i)*hexafib(i); dhexafib(i/2) = dhexafib(i/2) +kB(i)*hexafib(i); end

end

%destruction kinetics (fibres)

dmono = dmono +sel_red*2/i*kred*red*hexafib(i); dhexafib(i) = dhexafib(i) -sel_red*2/i*kred*red*hexafib(i); dred = dred -sel_red*2/i*kred*red*hexafib(i);

end

%destruction kinetics (trimers and tetramers)

dmono = dmono +(1-sel_red)*kred*red*tri +(1-sel_red)*kred*red*tetra; dtri = dtri -(1-sel_red)*kred*red*tri;

dtetra = dtetra -(1-sel_red)*kred*red*tetra; dred = dred -(1-sel_red)*kred*red*tri -(1-sel_red)*kred*red*tetra;

(18)

165 Appe

nd

ix

A. S

ou

rce

co

de

%calculating and displaying quench and redox fluxes

redflux=0; for i=2:maxlen redflux=redflux+kred*red*hexafib(i); end %{ if (mod(t,1e4)<10) disp(['t=',num2str(t)]); disp([' quench=',num2str(round(kquench*red*ox,5,'significant'))]); disp([' red=',num2str(round(redflux,5,'significant'))]); disp([' ox=',num2str(round(kox*ox*di+kox*ox*mono,5,'significant'))]); end %}

%storing of the differences in a column vector for the ode-solver

dy=zeros(n_ode,1); dy(1)=dmono; dy(2)=ddi; dy(3)=dtri; dy(4)=dtetra; dy(5)=dhexa; dy(6)=dred; dy(7)=dox; for i=2:maxlen dy(7+i-1)=dhexafib(i); end

A.2.3 Flux equation set

% ODEsim 1.0 % Fibres model % Elio Mattia % 25/3/2014

function dy=fluxFibres5_4(y,par)

%fetching maximum fibre length

maxlen = par(1);

%calculating number of ODEs

n_ode = 7+maxlen-1; %#ok<NASGU> %preallocating memory

kB=zeros(1,maxlen); hexafib=zeros(1,maxlen);

%fetching kinetic constants

kO2 = par(2); kred = par(3); kox = par(4); kquench = par(5); kflowox = par(6); kflowred = par(7); kX = par(8); kXF = par(9); kN = par(10); kS = par(11); kCE = par(12); for k=2:maxlen kB(k) = par(12+k-1); end fstart = par(12+maxlen); sel_red = par(12+maxlen+1); sel_ox = par(12+maxlen+2);

%no O2 oxidation for end flow computation

kO2 = 0; %fetching variables mono=y(1); di=y(2); tri=y(3); tetra=y(4); hexa=y(5); red=y(6); ox=y(7); for i=2:maxlen hexafib(i)=y(7+i-1); end

(19)

166 Appe

nd

ix

A. S

ou

rc

e c

od

e

%computation of fluxes %fred = 0; %initialized later

fox = kO2*mono +kO2*di +kO2*di +kox*ox*mono +kox*ox*di +kox*ox*di; fox34 = 2*(1-sel_ox)*kox*ox*di; fox6 = 2*sel_ox*kox*ox*di; fred34 = (1-sel_red)*kred*red*tri+(1-sel_red)*kred*red*tetra; fred6 = 0; fx = kX*tri*(mono+di) +kX*tetra*(mono+di); fxf = 0; %INACTIVE PATHWAY fn = 0; fs = 0; %INACTIVE PATHWAY fce = 0; fb = 0;

fflow = kflowred +kflowox;

fquench = kquench*red*ox +kquench*red*ox;

for i=2:maxlen if (i==2) %nucleation kinetics fn = fn +kN*tri*tri +kN*tetra*tetra; else %sequestration kinetics % INACTIVE PATHWAY

%catalyzed elongation kinetics

%catalyzed elongation flux

%only includes trimers and tetramers

fce = fce +kCE*1*tri*tri*2/(i-1)*hexafib(i-1) +kCE*1*tetra*tetra*2/(i-1)*hexafib(i-1); end

if (i>=4)

%fibre breakage kinetics

fb = fb +kB(i)*hexafib(i); end

%fibre destruction kinetics

fred6 = fred6 + sel_red*2/(i-1)*kred*red*hexafib(i);

end

fred = fred34+fred6;

%computation of total and relative fluxes

ftot= fred +fox +fx +fxf +fn +fs +fce +fb;

rfred = fred / ftot; rfox = fox / ftot; rfx = fx / ftot; rfxf = fxf / ftot; rfn = fn / ftot; rfs = fs / ftot; rfce = fce / ftot; rfb = fb / ftot;

rflowfflow = fflow / fflow; rflowfred = fred / fflow; rflowfox = fox / fflow; rflowfquench = fquench / fflow;

rfred34 = fred34 / (fred34+fred6); rfred6 = fred6 / (fred34+fred6); rfox34 = fox34 / (fox34+fox6); rfox6 = fox6 / (fox34+fox6);

rfredredox = (fred34+fred6) / (fred34+fred6+fox34+fox6); rfreplred = (fn+fce) / (fred34+fred6);

rfquenchred = fquench / (fquench+fred34+fred6);

%storing of the fluxes in a row vector

dy=[ftot, fflow, rfred, rfox, rfx, rfxf, rfn, rfs, rfce, rfb, rflowfflow, rflowfred, rflowfox, rflowfquench, rfred34, rfred6, rfox34, rfox6, rfredredox, rfreplred,rfquenchred];

% 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

A.2.4 3D plot toolset

% ODEsim 1.0 % Fibres model % Elio Mattia % 25/3/2014 %init clear all close all clc %Plot options

(20)

167 Appe

nd

ix

A. S

ou

rce

co

de

LineWidth = 2; FontSize = 22; %time tic

n_flux=14; %dimensions kflow=logspace(-10,-5,11); sel_red=linspace(0,1,11); sel_ox=linspace(0,1,11); %log dimensions lkflow=log10(kflow); lsel_red=sel_red; lsel_ox=sel_ox; %import fullStruct=load('kineticsFibres5_4_multiRun.mat'); flux=fullStruct.flux; sshexa=fullStruct.sshexa; %turnover processing replred=squeeze(flux(20,:,:,:)); quench=squeeze(flux(21,:,:,:)); turnover=(0.000008*8/24).*(1./replred).*(1./(1-quench)); for jkflow=1:length(kflow) turnover(jkflow,:,:)=turnover(jkflow,:,:)./kflow(jkflow); end nomturnover=0.*replred+1; for jkflow=1:length(kflow) nomturnover(jkflow,:,:)=nomturnover(jkflow,:,:).*(0.000008*8/24)./kflow(jkflow); end %flux selection psshexaconc=sshexa; pquench=squeeze(flux(21,:,:,:)); reltt=squeeze(bsxfun(@min,flux(15,:,:,:),flux(17,:,:,:))); ptt=reltt.*(1-pquench); relhexa=squeeze(bsxfun(@min,flux(16,:,:,:),flux(18,:,:,:))); phexa=relhexa.*(1-pquench); pcycle=squeeze(flux(20,:,:,:)); pturnover=turnover; %plot settings az=-37.5; el=+30.0;

axislimits=[lsel_red(1) lsel_red(length(lsel_red)) lkflow(1) lkflow(length(lkflow)) lsel_ox(1) lsel_ox(length(lsel_ox))];

i=0;

%ISOSURFACE COMPUTATION AND DISPLAY

i=i+1; figure(i) hFig = figure(i); set(hFig, 'Position', [50 100 1800 800]) plotrows=1; plotcolumns=2; %psshexaconc subplot(plotrows,plotcolumns,1)

set(gca, 'FontSize', FontSize, 'LineWidth', LineWidth); hold on

selectedflux=psshexaconc;

title('Steady state hexamer concentration') thresholds=[0.1 0.2 0.5 0.8 0.9];

%trace

p1 = patch(isosurface(lsel_red,lkflow,lsel_ox,selectedflux,thresholds(1))); isonormals(lsel_red,lkflow,lsel_ox,selectedflux,p1)

set(p1,'FaceColor','blue','EdgeColor','black');

set(p2,'FaceColor','green','EdgeColor','black');

set(p3,'FaceColor','yellow','EdgeColor','black');

set(p4,'FaceColor','red','EdgeColor','black');

set(p5,'FaceColor','black','EdgeColor','black'); daspect([1,1,1])

view(3); axis(axislimits);

(21)

168 Appe

nd

ix

A. S

ou

rc

e c

od

e

xlabel('x'); ylabel('log(kflow)'); zlabel('y'); axis square; view(az,el); grid on camlight lighting gouraud alpha(.5) %psshexaconc subplot(plotrows,plotcolumns,2)

selectedflux=psshexaconc;

title('Steady state hexamer concentration') thresholds=[0.1 0.2 0.5 0.8 0.9];

%trace

set(p5,'FaceColor','black','EdgeColor','black'); daspect([1,1,1]) view(3); axis(axislimits); xlabel('x'); ylabel('log(kflow)'); zlabel('y'); camlight lighting gouraud axis square; axis ij; view(az,el); grid on alpha(.5)

i=i+1; figure(i) hFig = figure(i); set(hFig, 'Position', [50 100 1800 800]) plotrows=1; plotcolumns=2; %pquench subplot(plotrows,plotcolumns,1)

selectedflux=pquench;

title('Steady state quench flux') thresholds=[0.1 0.2 0.5 0.8 0.9];

%trace

set(p5,'FaceColor','black','EdgeColor','black'); daspect([1,1,1]) view(3); axis(axislimits); xlabel('x'); ylabel('log(kflow)'); zlabel('y'); camlight

(22)

169 Appe

nd

ix

A. S

ou

rce

co

de

lighting gouraud axis square; view(az,el); grid on alpha(.5) %pquench subplot(plotrows,plotcolumns,2)

selectedflux=pquench;

title('Steady state quench flux') thresholds=[0.1 0.2 0.5 0.8 0.9];

%trace

set(p5,'FaceColor','black','EdgeColor','black'); daspect([1,1,1]) view(3); axis(axislimits); xlabel('x'); ylabel('log(kflow)'); zlabel('y'); camlight lighting gouraud axis square; axis ij; view(az,el); grid on alpha(.5)

i=i+1; figure(i) hFig = figure(i); set(hFig, 'Position', [50 100 1800 800]) plotrows=1; plotcolumns=2; %ptt subplot(plotrows,plotcolumns,1)

selectedflux=ptt;

title('Steady state trimers/tetramers recycling flux') thresholds=[0.1 0.2 0.5 0.8 0.9];

%trace

set(p5,'FaceColor','black','EdgeColor','black'); daspect([1,1,1]) view(3); axis(axislimits); xlabel('x'); ylabel('log(kflow)'); zlabel('y'); camlight lighting gouraud axis square; view(az,el); grid on alpha(.5)