1PySCeS - Python Simulator for Cellular Systems (http://pysces.sourceforge.net)
2"""
3
4Copyright (C) 2004-2020 B.G. Olivier, J.M. Rohwer, J.-H.S Hofmeyr all rights reserved,
5
6Brett G. Olivier (bgoli@users.sourceforge.net)
7Triple-J Group for Molecular Cell Physiology
8Stellenbosch University, South Africa.
9
10Permission to use, modify, and distribute this software is given under the
11terms of the PySceS (BSD style) license. See LICENSE.txt that came with
12this distribution for specifics.
13
14NO WARRANTY IS EXPRESSED OR IMPLIED.  USE AT YOUR OWN RISK.
15Brett G. Olivier
16"""
17from __future__ import division, print_function
18from __future__ import absolute_import
19from __future__ import unicode_literals
20
21from pysces.version import __version__
22__doc__ = '''SBML reading/writing module - now replaced by PySCeS Core2'''
23
24import os,sys
25from time import sleep, strftime
26from getpass import getuser
27try:
28    input = raw_input  # Py2 compatibility
29except NameError:
30    pass
31
32if sys.platform == 'win32':
33    try:
34        import pysces.libsbml.libsbmlinit as SBML
35    except Exception as e:
36        print('Windows sbml load error',e)
37else:
38    try:
39        import libsbml as SBML
40    except Exception as e:
41        print('Posix sbml load error',e)
42
43
44class PyscesSBML:
45    """The PySCeS interface to libSBML and SBML utilities"""
46
47    mode_number_format = '%2.4e'
48    sbml_level = 2
49    _debug = 0
50    SBML = SBML
51
52    def SBML_buildBasicModel(self,mod,filename,slvl=2,dir=None,substance=(1,0),volume=(1,0),time=(1,0),arules=None,notes=None):
53        """
54        SBML_buildBasicModel(mod,filename,slvl=2,dir=None)
55
56        Create a basic SBML model.
57        Arguments:
58        =========
59        mod: an active PySCeS model object
60        filename: the output SBML file name
61        slvl [default=2]: the SBML level that should be used
62        dir [default=None]: the output directory
63        substance [default=(1,0)]: the model substance unit - SBML default is "mole"
64        volume [default=(1,0)]: the model volume unit - SBML default is "litre"
65        time [default=(1,0)]: the model time unit - SBML default is "second"
66
67        """
68        self.SBML_createModel(mod,filename,slvl,dir)
69        self.SBML_setCompartment()
70        self.SBML_setNotes(txt=notes)
71        self.SBML_setUnits(substance=substance,volume=volume,time=time)
72        if arules != None: self.SBML_setAssignmentRules(arules)
73        self.SBML_setSpecies()
74        self.SBML_setReactions()
75        self.SBML_setModel()
76
77
78    def write(self,mod,filename,slvl=2,dir=None,substance=(1,0),volume=(1,0),time=(1,0),arules=None,notes=None):
79        """
80        write(mod,filename,slvl=2,dir=None)
81
82        Write a PySCeS model as an SBML file.
83
84        Arguments:
85        =========
86        mod: an active PySCeS model object
87        filename: the output SBML file name
88        slvl [default=2]: the SBML level that should be used
89        dir [default=None]: the output directory
90        substance [default=(1,0)]: the model substance unit - SBML default is "mole"
91        volume [default=(1,0)]: the model volume unit - SBML default is "litre"
92        time [default=(1,0)]: the model time unit - SBML default is "second"
93
94        """
95        self.SBML_buildBasicModel(mod,filename,slvl,dir,substance,volume,time,arules,notes)
96        self.SBML_writeFile()
97
98
99    def getSBML_document(self,mod,substance=(1,0),volume=(1,0),time=(1,0),arules=None,notes=None):
100        """
101        Returns an SBML document object
102
103        Arguments:
104        =========
105        mod: an active PySCeS model object
106        substance [default=(1,0)]: the model substance unit - SBML default is "mole"
107        volume [default=(1,0)]: the model volume unit - SBML default is "litre"
108        time [default=(1,0)]: the model time unit - SBML default is "second"
109
110        """
111
112        filename = 'tempXML'
113        slvl = 2
114        dir = None
115        self.SBML_buildBasicModel(mod,filename,slvl,dir,substance,volume,time,arules,notes)
116        return self.sbml_document
117
118    def getSBML_string(self,mod,substance=(1,0),volume=(1,0),time=(1,0),arules=None,notes=None):
119        """
120        Returns an SBML file as a string
121
122        Arguments:
123        =========
124        mod: an active PySCeS model object
125        substance [default=(1,0)]: the model substance unit - SBML default is "mole"
126        volume [default=(1,0)]: the model volume unit - SBML default is "litre"
127        time [default=(1,0)]: the model time unit - SBML default is "second"
128        """
129
130        filename = 'tempXML'
131        slvl = 2
132        dir = None
133        self.SBML_buildBasicModel(mod,filename,slvl,dir,substance,volume,time,arules,notes)
134        return self.sbml_document.toSBML()
135
136    def __cleanString__(self,s):
137        s = s.lstrip()
138        s = s.rstrip()
139        return s
140
141    def parseForcingFunctions(self):
142        self.__forcing_function_dic__ = {}
143        ff = self.model_obj._Function_forced.split('\n')
144        for f in ff:
145            if f != '':
146                f =  f.split('=')
147                f[0] = f[0].replace('self.','')
148                f[1] = f[1].replace('self.','')
149                self.__forcing_function_dic__.setdefault(self.__cleanString__(f[0]), self.__cleanString__(f[1]))
150
151    def getAssignmentRules(self):
152        self.parseForcingFunctions()
153        out = []
154        for key in list(self.__forcing_function_dic__.keys()):
155            out.append((key, self.__forcing_function_dic__[key]))
156        return out
157
158
159
160    def SBML_createModel(self,mod,filename,slvl=2,dir=None):
161        """
162        SBML_createModel(mod,filename,slvl=2,dir=None)
163
164        Set up an SBML document and extract the model NetworkDict
165
166        Arguments:
167        =========
168        mod: a PySCeS model object
169        filename: the output filename
170        slvl [default=2]: SBML level required
171        dir [default=None]: output directory
172
173        """
174        if self._debug: print('SBML_createModel')
175        self.model_obj = mod
176        self.__nDict__= self.model_obj.__nDict__
177        self.model_filename = filename
178        if dir == None:
179            self.model_dir = self.model_obj.ModelOutput
180        else:
181            self.model_dir = dir
182
183        self.sbml_level = slvl
184        self.sbml_model = self.SBML.Model()
185        self.sbml_model.setName(self.model_obj.ModelFile[:-4])
186        self.sbml_document = self.SBML.SBMLDocument()
187        # new stuff
188        self.global_parameters = []
189        self.model_compartment_name = None
190
191        # create initdict
192        self.__InitStrings__ = [s.replace('self.','') for s in self.model_obj._PysMod__InitStrings]
193        self.__InitDict__ = {}
194        for ii in self.__InitStrings__:
195            l,r = ii.split('=')
196            self.__InitDict__.setdefault(self.__cleanString__(l), float(self.__cleanString__(r)))
197
198        # create forcing function dic
199        try:
200            self.parseForcingFunctions()
201        except:
202            print("No pre-defined forcing functions")
203
204
205        if self.sbml_level == 1:
206            if sys.platform == 'win32':
207                print('Due to a bug in self.SBML for Windows writing a lvl 1 file will crash your session writing lvl 2 instead ... sorry')
208                self.sbml_document.setLevel(2)
209            else:
210                self.sbml_document.setLevel(self.sbml_level)
211        else:
212            self.sbml_document.setLevel(2)
213
214    def SBML_setCompartment(self,name=None,vol=1):
215        """
216        SBML_setCompartment(name=None,vol=1)
217
218        Initialise SBML compartments (note PySCeS currently utilises a single compartment)
219
220        Arguments:
221        =========
222        name [default=None]: the compartment name, default is compartment1
223        vol [default=1]: the compartment volume
224
225        """
226        if self._debug: print('SBML_setCompartment')
227        comp_def = self.sbml_model.createCompartment()
228        if not name:
229            self.model_compartment_name = 'Cell'
230        else:
231            self.model_compartment_name = name
232            for char in [' ','.','-','*','?','!','\t','\n']:
233                self.model_compartment_name = self.model_compartment_name.replace(char,'_')
234            self.model_compartment_name = name
235        comp_def.setId(self.model_compartment_name)
236        comp_def.setVolume(vol)
237
238    def SBML_setNotes(self,txt=None):
239        notes = '<body xmlns="http://www.w3.org/1999/xhtml">'
240        if txt != None:
241            notes += '<span style="font-family: Courier New,Courier,monospace;">'
242            notes += txt
243            notes += '</span>'
244        notes += '</body>'
245        self.sbml_model.setNotes(notes)
246
247
248    def SBML_setUnits(self, **kwargs):
249        """
250        SBML_setUnits(substance=(1,0), volume=(1,0), time=(1,0))
251
252        Set the SBML default units note that the input here is the factor and index multiplying
253        the SBML default, so for example the default substance (1,0) is (1*10**0)*mole So e.g.
254        if you were specifing default units of millimoles and seconds you would
255        set substance=(1,-3) and time=(60,0) i.e. (1*10**-3)*mole and (60*10**0)*seconds
256
257        Arguments:
258        =========
259        substance [default=(1,0)]: the model substance unit - SBML default is "mole"
260        volume [default=(1,0)]: the model volume unit - SBML default is "litre"
261        time [default=(1,0)]: the model time unit - SBML default is "second"
262
263        """
264        for un in list(kwargs.keys()):
265            vdef = self.sbml_model.createUnitDefinition()
266            vdef.setId(un)
267            vu = self.sbml_model.createUnit()
268            if un == 'substance': vu.setKind(self.SBML.UnitKind_forName('mole'))
269            elif un == 'volume': vu.setKind(self.SBML.UnitKind_forName('litre'))
270            elif un == 'time': vu.setKind(self.SBML.UnitKind_forName('second'))
271            vu.setMultiplier(kwargs[un][0])
272            vu.setScale(kwargs[un][1])
273            vu.setOffset(0)
274
275
276    def SBML_setSpecies(self):
277        """
278        SBML_setSpecies()
279
280        Initialise and add species information to the SBML model
281
282        Arguments:
283        None
284
285        """
286        if self._debug: print('SBML_setSpecies')
287        reagList = self.model_obj.__species__ + self.model_obj.__fixed_species__
288        for reagent in range(len(reagList)):
289            s = self.sbml_model.createSpecies()
290            s.setId(reagList[reagent])
291            s.setName(reagList[reagent])
292            s.setCompartment(self.model_compartment_name)
293            if reagList[reagent] in self.model_obj.__fixed_species__:
294                s.setBoundaryCondition(True)
295                s.setConstant(True)
296            else:
297                s.setBoundaryCondition(False)
298
299            if reagent < len(self.model_obj.__species__ ):
300                reagName = reagList[reagent] + '_init'
301            else:
302                reagName = reagList[reagent]
303
304            if self.sbml_level == 1:
305                s.setInitialAmount(getattr(self.model_obj,reagName))
306            else:
307                s.setInitialConcentration(getattr(self.model_obj,reagName))
308
309    def SBML_setAssignmentRules(self, rules=[]):
310
311        for rule in rules:
312            print(rule)
313            self.global_parameters.append(rule[0])
314            p = self.sbml_model.createParameter()
315            p.setId(rule[0])
316            p.setValue(getattr(self.model_obj,rule[0]))
317            p.setConstant(False)
318            r = self.sbml_model.createAssignmentRule()
319            r.setVariable(rule[0])
320            r.setFormula(rule[1])
321            r.setMathFromFormula()
322
323
324    def SBML_setReactions(self):
325        """
326        SBML_setReactions()
327
328        Add kinetic rate laws to the SBMl model
329
330        Arguments:
331        None
332
333        """
334        if self._debug: print('SBML_setReactions')
335        #TotSpecies = list(self.model_obj._PysMod__FixedReagents)+list(self.model_obj._PysMod__VarReagents)
336        reaction_params = []
337        for rxn in self.model_obj._PysMod__ReactionIDs:
338            print('Adding reaction:', rxn)
339            i = self.sbml_model.createReaction()
340            i.setId(rxn)
341            ndr = self.model_network_dict[rxn]
342            for reagent in ndr['Reagents']:
343                stoich = ndr['Reagents'][reagent]
344                species = self.SBML.SpeciesReference(reagent.replace('self.',''),abs(stoich))
345                if stoich < 0:
346                    i.addReactant(species)
347                elif stoich > 0:
348                    i.addProduct(species)
349                elif stoich == 0:
350                    i.addModifier(species)
351            # add a volume to convert rate equation to kinetic law
352            kineticLaw = ndr['RateEq'].replace('self.','')
353            kineticLaw = kineticLaw.replace('scipy.','')
354            if self.model_compartment_name not in self.model_obj.parameters:
355                kineticLaw = self.model_compartment_name + ' * ('+ kineticLaw +')'
356            else:
357                kineticLaw = kineticLaw
358            kineticLaw = self.SBML.KineticLaw(kineticLaw)
359
360            # local parameters retired in favour of globals
361            ##  for parameter in ndr['Params']:
362                ##  p = parameter.replace('self.','')
363                ##  if p not in self.model_obj.__fixed_species__ and p not in self.global_parameters:
364                    ##  try:
365                        ##  kineticLaw.addParameter(self.SBML.Parameter(p, getattr(self.model_obj,p)))
366                        ##  reaction_params.append(p)
367                    ##  except AttributeError,err :
368                        ##  print '\n', err
369                        ##  print "Parameter set error ... are there forcing functions??"
370                        ##  sleep(0.5)
371            i.setKineticLaw(kineticLaw)
372            if ndr['Type'] == 'Rever':
373                rev = True
374            else:
375                rev = False
376            i.setReversible(rev)
377
378            # Add modifiers to reaction - brett 20050607
379            for reac in self.model_obj.__modifiers__:
380                if reac[0] == rxn:
381                    for x in reac[1]:
382                        print(' ' + reac[0] +' has modifier: ' + x)
383                        self.sbml_model.createModifier().setSpecies(x)
384
385        # add extra parameter initialised but not in reactions
386        # we have to do this in case the assignment rules are added after we build the model
387        hack = list(self.__forcing_function_dic__.keys())
388
389        not_xparams = self.global_parameters + reaction_params+\
390                     list(self.model_obj.species)+\
391                     list(self.model_obj.fixed_species) + [self.model_compartment_name] +hack
392
393
394        for k in list(self.__InitDict__.keys()):
395            if k not in not_xparams:
396                print('Adding parameter:', k)
397                self.global_parameters.append(k)
398                p = self.sbml_model.createParameter()
399                p.setId(k)
400                p.setValue(getattr(self.model_obj, k))
401
402
403
404
405    def SBML_setModel(self):
406        """
407        SBML_setModel()
408
409        Add the SBML model to the predefined SBML document
410
411        Arguments:
412        None
413
414        """
415        if self._debug: print('SBML_setModel')
416        self.sbml_document.setModel(self.sbml_model)
417
418    def SBML_writeFile(self):
419        """
420        SBML_writeFile()
421
422        Write the SBML document to predefined output file
423
424        Arguments:
425        None
426
427        """
428
429        self.SBML.writeSBML(self.sbml_document,'pysces_sbml_tmp.xml')
430        Fin = open('pysces_sbml_tmp.xml','r')
431        Fout = open(os.path.join(self.model_dir,self.model_filename+'.xml'),'w')
432        cntr = 0
433        try:
434            UseR = getuser()
435        except:
436            UseR = ''
437        for line in Fin:
438            if cntr == 1:
439                Fout.write('<!-- Created with PySCeS ('+ __version__ + ') on ' + strftime("%a, %d %b %Y %H:%M:%S") + ' by '+UseR+' -->\n'+line)
440            else:
441                Fout.write(line)
442            cntr += 1
443        Fout.close()
444        Fin.close()
445
446        os.remove('pysces_sbml_tmp.xml')
447
448
449    def convert2psc(self,filename,dir=None,dirOut=None):
450        """
451        convert2psc(filename,dir=None,dirOut=None)
452
453        Convert an SBML file into a PySCeS input file
454
455        Arguments:
456        =========
457        filename: the SBML source file
458        dir [default=None]: specify the SBMl file directory
459        dirOut [default=None]: the PSC file output directory
460
461        """
462        if dir == None:
463            dir = os.getcwd()
464        File = os.path.join(dir,filename)
465        assert os.path.exists(File), "Invalid path"
466
467        self.model_filename = filename
468
469        r = self.SBML.SBMLReader()
470        d = r.readSBML(File)
471        m = d.getModel()
472
473        def getName(i):
474            if d.getLevel() == 1:
475                return i.getName()
476            else:
477                return i.getId()
478
479        reactions = m.getListOfReactions()
480        ReactionIDs = []
481        for i in reactions:
482            ReactionIDs.append(getName(i))
483
484        init_fixed = []
485        init_var = []
486        init_par = []
487        parameters = []
488
489        for i in m.getListOfSpecies():
490            parName = getName(i)
491        # if a species is a BoundaryCondition or constant it becomes fixed - brett 20050111
492            if i.getBoundaryCondition() or i.getConstant():
493                if i.getConstant() and not i.getBoundaryCondition():
494                    print(parName, ' is set as constant, assuming: BoundaryCondition = True')
495                init_fixed.append((parName,i.getInitialConcentration()))
496            else:
497                init_var.append((parName,i.getInitialConcentration()))
498
499        NetworkDict = dict([(i,dict.fromkeys(['Params',
500                              'RateEq',
501                              'Reagents',
502                              'Type'])) for i in ReactionIDs])
503
504
505        for i in reactions:
506            rDict = NetworkDict[getName(i)]
507            j = i.getKineticLaw()
508            par = []
509            try:
510                for k in j.getListOfParameters():
511                    par.append(getName(k))
512                    init_par.append((getName(k),k.getValue()))
513                    parameters.append(getName(k))
514                rDict['Params'] = par
515                rDict['RateEq'] = j.getFormula()
516                if d.getLevel() == 1:
517                    rDict['RateEq'] = rDict['RateEq'].replace(' ', '' )
518                    rDict['RateEq'] = rDict['RateEq'].replace('^', '**')
519            except Exception as err:
520                rDict['Params'] = []
521                rDict['RateEq'] = ''
522                print(err)
523
524            Substrates = []
525            Products = []
526
527            for k in i.getListOfReactants():
528                species = k.getSpecies()
529                stoich = -k.getStoichiometry()
530                Substrates.append((species,stoich))
531
532            for k in i.getListOfProducts():
533                species = k.getSpecies()
534                stoich = k.getStoichiometry()
535                Products.append((species,stoich))
536
537            # this is to eliminate zero stoichiometries {0}xyz
538            badList = []
539            for sub in Substrates:
540                if sub[1] == 0:
541                    badList.append(sub)
542            for bad in badList:
543                Substrates.pop(Substrates.index(bad))
544
545            badList = []
546            for prod in Products:
547                if prod[1] == 0:
548                    badList.append(prod)
549            for bad in badList:
550                Products.pop(Products.index(bad))
551
552            # add source/sink pools to nasty substrate/productless reactions - brett 20050908
553            if len(Substrates) == 0:
554                Substrates.append(('$pool',-1.0))
555            if len(Products) == 0:
556                Products.append(('$pool',1.0))
557
558            # print Substrates
559            # print Products
560
561            rDict['Reagents'] = dict(Substrates+Products)
562            if i.getReversible() == True:
563                t = 'Rever'
564            else:
565                t = 'Irrev'
566            rDict['Type'] = t
567            NetworkDict[getName(i)].update(rDict)
568
569        # Add extra model parameters not defined in reactions (apparently)
570        if len(m.getListOfParameters()) > 0:
571            for x in m.getListOfParameters():
572                if getName(x) not in parameters:
573                    #print getName(x)
574                    init_par.append((getName(x),x.getValue()))
575
576        if dirOut == None:
577            self.model_filename = os.path.join(os.getcwd(),self.model_filename)
578        else:
579            self.model_filename = os.path.join(dirOut,self.model_filename)
580
581        #print 'init_par'
582        #print init_par
583        #print 'init_var'
584        #print init_var
585        #print 'init_fixed'
586        #print init_fixed
587
588        # sometimes things just work lekker (replaced all the old showS^&t) - brett 20050913
589        outFile = open(self.model_filename+'.psc','w')
590        self.PSC_writeHeader(outFile)
591        self.PSC_writeFixedSpeciesList(outFile,init_fixed)
592        self.PSC_writeRateEquations(outFile,NetworkDict,number_format='%2.3f')
593        self.PSC_writeSpecies(outFile,init_var)
594        self.PSC_writeFixedSpecies(outFile,init_fixed)
595        self.PSC_writeParameters(outFile,init_par)
596        outFile.close()
597
598        # Initialise compartment volumes as a parameter - brett 20050908
599        compartmentList = []
600        for comp in m.getListOfCompartments():
601            #print comp
602            compartmentList.append((getName(comp),comp.getVolume()))
603
604        if len(compartmentList) > 1:
605            print('\nINFO: PySCeS models are assumed to have a single compartment')
606
607        if len(compartmentList) > 0:
608            F = open(self.model_filename+'.psc','a')
609            F.write('\n## Initialise compartment volumes')
610            for comp in compartmentList:
611                F.write('\n' + comp[0] + ' = ' + str(comp[1]))
612                ##  parameters.append(x[0])
613                ##  init_par.append(x)
614            F.write('\n')
615            F.close()
616
617
618        # Add assignment rules as forcing functions - brett 20050908
619        pscRules = []
620        for rule in m.getListOfRules():
621            pscRules.append((rule.getVariable(),rule.getFormula()))
622
623        if len(pscRules) > 0:
624            F = open(self.model_filename+'.psc','a')
625            F.write('\n## Assignment rules translated to forcing functions\n')
626            for rule in pscRules:
627                rule0 = 'self.' + rule[0]
628                rule1l = rule[1].split()
629                for word in range(len(rule1l)):
630                    if rule1l[word].isalnum():
631                        if rule1l[word] not in ['1','2','3','4','5','6','7','8','9']:
632                            rule1l[word] = 'self.'+rule1l[word]
633                F.write('!F '+ rule0 + ' = ')
634                for word in rule1l:
635                    F.write(word + ' ')
636                F.write('\n')
637            F.write('\n')
638            F.close()
639
640        if len(m.getNotes()) > 0:
641            F = open(self.model_filename+'.psc','a')
642            F.write('\n## Model notes' + m.getNotes().replace('\n','\n# ')+'\n\n')
643            F.close()
644
645
646    def PSC_writeHeader(self,File):
647        """
648        PSC_writeHeader(File)
649
650        Write a PSC file header to an open file object
651
652        Arguments:
653        =========
654        File: a writable open text file object
655
656        """
657        try:
658            UseR = getuser()
659        except:
660            UseR = ''
661
662        header = ''
663        #header +=  '############################################################\n'
664        header += '# PySCeS (' + __version__ + ') model input file\n'
665        header +=  '# PySCeS can be found at http://pysces.sourceforge.net/\n'
666        #header += '###########################################################\n\n'
667        header += '# Original input file: ' + File.name.split('\\')[-1][:-4] + '\n'
668        header += '# This file generated: ' + strftime("%a, %d %b %Y %H:%M:%S") + ' by '+UseR+'\n\n'
669        File.write(header)
670        File.write('\n')
671
672
673    def PSC_writeSpecies(self,File,species):
674        """
675        PSC_writeSpecies(File,species)
676
677        Write out model species initiaiisations to file
678
679        Arguments:
680        =========
681        File: a writable open file object
682        species: a list of (species.value) pairs
683
684        """
685        out_list = []
686
687        out_list.append('\n## Variable species initial values\n')
688        for x in range(len(species)):
689            out_list.append(species[x][0] + ' = ' +  self.mode_number_format % species[x][1] + '\n')
690        for x in out_list:
691            File.write(x)
692        File.write('\n')
693
694
695    def PSC_writeFixedSpeciesList(self,File,fixed_species):
696        """
697        PSC_writeFixedSpeciesList(File,fixed_species)
698
699        Write fixed species declaration to a PSC file
700
701        Arguments:
702        =========
703        File: open, writable file object
704        fixed_species: a list of (species,value) pairs
705
706        """
707        File.write('## Fixed species\n')
708        if len(fixed_species) == 0:
709            File.write('#  <none>')
710        else:
711            File.write('FIX: ')
712            for x in fixed_species:
713                File.write(x[0] + ' ')
714        File.write('\n\n')
715
716
717    def PSC_writeFixedSpecies(self,File,fixed_species):
718        """
719        PSC_writeFixedSpecies(File,fixed_species)
720
721        Write fixed species initialisations to a PSC file
722
723        Arguments:
724        =========
725        File: open, writable file object
726        fixed_species: a list of (species,value) pairs
727
728        """
729        out_list = []
730
731        out_list.append('\n## Fixed species\n')
732        for x in range(len(fixed_species)):
733            out_list.append(fixed_species[x][0] + ' = ' +  self.mode_number_format % fixed_species[x][1] + '\n')
734        for x in out_list:
735            File.write(x)
736        File.write('\n')
737
738
739    def PSC_writeParameters(self,File,parameters):
740        """
741        PSC_writeParameters(File,parameters)
742
743        Write mode parameter initialisations to a PSC file
744
745        Arguments:
746        =========
747        File: open, writable file object
748        parameters: a list of (parameter,value) pairs
749
750        """
751        out_list = []
752
753        out_list.append('\n## Parameters\n')
754        for x in range(len(parameters)):
755            out_list.append(parameters[x][0] + ' = ' +  self.mode_number_format % parameters[x][1] + '\n')
756        for x in out_list:
757            File.write(x)
758        File.write('\n')
759
760
761    def PSC_writeRateEquations(self,File,NetworkDict,number_format='%2.3f'):
762        """
763        PSC_writeRateEquations(File,NetworkDict,number_format='%2.3f')
764
765        Write model rate equations to a PSC file
766
767        Arguments:
768        =========
769        File: open, writable file object
770        NetworkDict: a PySCeS network dictionary
771        number_format [default='%2.3f']: number formatting to use in rate laws
772
773        """
774        out_list = []
775
776        out_list.append('\n## Reaction stoichiometry and rate equations\n')
777        for key in NetworkDict:
778            out_list.append(key + ':\n')
779            reagL = []
780            reagR = []
781            for reagent in NetworkDict[key]['Reagents']:
782                if NetworkDict[key]['Reagents'][reagent] > 0:
783                    if NetworkDict[key]['Reagents'][reagent] == 1.0:
784                        reagR.append(reagent.replace('self.',''))
785                    else:
786                        reagR.append('{' + number_format % abs(NetworkDict[key]['Reagents'][reagent]) + '}' + reagent.replace('self.',''))
787                elif NetworkDict[key]['Reagents'][reagent] < 0:
788                    if NetworkDict[key]['Reagents'][reagent] == -1.0:
789                        reagL.append(reagent.replace('self.',''))
790                    else:
791                        reagL.append('{' + number_format % abs(NetworkDict[key]['Reagents'][reagent]) + '}' + reagent.replace('self.',''))
792                elif NetworkDict[key]['Reagents'][reagent] == 0:
793                    #reagL.append(reagent.replace('self.',''))
794                    print(NetworkDict[key]['Reagents'])
795                    input('WTF: please contact developers')
796
797            if len(reagL) == 0:
798                print('Zero pool substrate', File.name)
799                reagL.append('$pool')
800            if len(reagR) == 0:
801                print('Zero pool product', File.name)
802                reagR.append('$pool')
803
804            substring = ''
805            count = 0
806            for x in reagL:
807                if count != 0:
808                    substring += ' + '
809                substring += x.replace(' ','')
810                count += 1
811            prodstring = ''
812            count = 0
813            for x in reagR:
814                if count != 0:
815                    prodstring += ' + '
816                prodstring += x.replace(' ','')
817                count += 1
818
819            if NetworkDict[key]['Type'] == 'Rever':
820                symbol = ' = '
821            else:
822                symbol = ' > '
823            out_list.append('\t' + substring + symbol + prodstring + '\n')
824            out_list.append('\t' + NetworkDict[key]['RateEq'].replace('self.','') + '\n\n')
825        for x in out_list:
826            File.write(x)
827