diff --git a/generateFluxDiagram.py b/generateFluxDiagram.py
index 28e6d4047b..be303863a9 100644
--- a/generateFluxDiagram.py
+++ b/generateFluxDiagram.py
@@ -47,7 +47,7 @@
 
 ################################################################################
 
-def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, outputDirectory, speciesDirectory=None, settings = None):
+def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, outputDirectory, centralSpecies = None, speciesDirectory=None, settings = None):
     """
     For a given `reactionModel` and simulation results stored as arrays of
     `times`, species `concentrations`, and `reactionRates`, generate a series
@@ -55,7 +55,7 @@ def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, out
     a movie. The individual frames and the final movie are saved on disk at
     `outputDirectory.`
     """
-    
+    global maximumNodeCount, maximumEdgeCount, timeStep, concentrationTolerance, speciesRateTolerance
     # Allow user defined settings for flux diagram generation if given
     if settings:
         maximumNodeCount = settings['maximumNodeCount']       
@@ -70,6 +70,13 @@ def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, out
     reactionList = reactionModel.core.reactions[:]
     numReactions = len(reactionList)
     
+    #search index of central species:
+    if centralSpecies is not None:
+        for i, species in enumerate(speciesList):
+            if species.label == centralSpecies:
+                centralSpeciesIndex = i
+                break 
+    
     # Compute the rates between each pair of species (big matrix warning!)
     speciesRates = numpy.zeros((len(times),numSpecies,numSpecies), numpy.float64)
     for index, reaction in enumerate(reactionList):
@@ -92,22 +99,40 @@ def generateFluxDiagram(reactionModel, times, concentrations, reactionRates, out
     
     # Determine the nodes and edges to keep
     nodes = []; edges = []
-    for i in range(numSpecies*numSpecies):
-        productIndex, reactantIndex = divmod(speciesIndex[-i-1], numSpecies)
-        if reactantIndex > productIndex:
-            # Both reactant -> product and product -> reactant are in this list,
-            # so only keep one of them
-            continue
-        if maxSpeciesRates[reactantIndex, productIndex] == 0:
-            break
-        if reactantIndex not in nodes and len(nodes) < maximumNodeCount: nodes.append(reactantIndex)
-        if productIndex not in nodes and len(nodes) < maximumNodeCount: nodes.append(productIndex)
-        if len(nodes) > maximumNodeCount: 
-            break
-        edges.append([reactantIndex, productIndex])
-        if len(edges) >= maximumEdgeCount:
-            break
-    
+    if centralSpecies is None:
+        for i in range(numSpecies*numSpecies):
+            productIndex, reactantIndex = divmod(speciesIndex[-i-1], numSpecies)
+            if reactantIndex > productIndex:
+                # Both reactant -> product and product -> reactant are in this list,
+                # so only keep one of them
+                continue
+            if maxSpeciesRates[reactantIndex, productIndex] == 0:
+                break
+            if reactantIndex not in nodes and len(nodes) < maximumNodeCount: nodes.append(reactantIndex)
+            if productIndex not in nodes and len(nodes) < maximumNodeCount: nodes.append(productIndex)
+            if len(nodes) > maximumNodeCount: 
+                break
+            edges.append([reactantIndex, productIndex])
+            if len(edges) >= maximumEdgeCount:
+                break
+    else:
+        nodes.append(centralSpeciesIndex)
+        for index, reaction in enumerate(reactionList):
+            for reactant, product in reaction.pairs:
+                reactantIndex = speciesList.index(reactant)
+                productIndex = speciesList.index(product)
+                if maxSpeciesRates[reactantIndex, productIndex] == 0:
+                    break
+                if len(nodes) > maximumNodeCount or len(edges) >= maximumEdgeCount: 
+                    break
+                if reactantIndex == centralSpeciesIndex: 
+                    if productIndex not in nodes:
+                        nodes.append(productIndex)
+                        edges.append([reactantIndex, productIndex])
+                if productIndex == centralSpeciesIndex: 
+                    if reactantIndex not in nodes:
+                        nodes.append(reactantIndex)
+                        edges.append([reactantIndex, productIndex])
     # Create the master graph
     # First we're going to generate the coordinates for all of the nodes; for
     # this we use the thickest pen widths for all nodes and edges 
@@ -224,7 +249,7 @@ def simulate(reactionModel, reactionSystem, settings = None):
     Generate and return a set of core and edge species and reaction fluxes
     by simulating the given `reactionSystem` using the given `reactionModel`.
     """
-    
+    global maximumNodeCount, maximumEdgeCount, timeStep, concentrationTolerance, speciesRateTolerance
     # Allow user defined settings for flux diagram generation if given
     if settings:
         maximumNodeCount = settings['maximumNodeCount']       
@@ -487,7 +512,7 @@ def loadRMGPyJob(inputFile, chemkinFile=None, speciesDict=None):
 
 ################################################################################
 
-def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = False, settings = None, chemkinOutput = ''):
+def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = False, settings = None, chemkinOutput = '', centralSpecies = None):
     """
     Generates the flux diagram based on a condition 'inputFile', chemkin.inp chemkinFile,
     a speciesDict txt file, plus an optional chemkinOutput file.
@@ -510,7 +535,7 @@ def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = Fals
         time, coreSpeciesConcentrations, coreReactionRates, edgeReactionRates = loadChemkinOutput(chemkinOutput, rmg.reactionModel)
 
         print 'Generating flux diagram for chemkin output...'
-        generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, os.path.join(savePath, '1'), speciesPath, settings)
+        generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, os.path.join(savePath, '1'), centralSpecies, speciesPath, settings)
 
     else:
         # Generate a flux diagram video for each reaction system
@@ -532,9 +557,8 @@ def createFluxDiagram(savePath, inputFile, chemkinFile, speciesDict, java = Fals
             time, coreSpeciesConcentrations, coreReactionRates, edgeReactionRates = simulate(rmg.reactionModel, reactionSystem, settings)
 
             print 'Generating flux diagram for reaction system {0:d}...'.format(index+1)
-            generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, 
-                                os.path.join(savePath, '{0:d}'.format(index+1)), speciesPath, settings)
-
+            generateFluxDiagram(rmg.reactionModel, time, coreSpeciesConcentrations, coreReactionRates, os.path.join(savePath, '{0:d}'.format(index+1)), 
+                                centralSpecies, speciesPath, settings)
 ################################################################################
 
 if __name__ == '__main__':