diff --git a/.gitignore b/.gitignore
index 9b703336..0fd85fa8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -173,3 +173,5 @@ arm64
 /pyneuroml/xppaut/report.txt
 /examples/test_data/xppaut/fhn_2.ode
 /examples/test_data/xppaut/fhn_brain.py
+/examples/test_data/xppaut/nca_2.ode
+/examples/test_data/xppaut/wc_2.ode
diff --git a/pyneuroml/analysis/LEMS_Test_TEMPLATE.xml b/pyneuroml/analysis/LEMS_Test_TEMPLATE.xml
index b131a1c7..5609ec96 100644
--- a/pyneuroml/analysis/LEMS_Test_TEMPLATE.xml
+++ b/pyneuroml/analysis/LEMS_Test_TEMPLATE.xml
@@ -139,10 +139,10 @@
 #end## 
 
 #foreach ($g in $gates)##
-        <OutputFile id="${g}_inf_of1" fileName="${channel}.${g}.inf.lems${dat_suffix}.dat">
+        <OutputFile id="${g}_inf_of1" fileName="${channel}.${g.replace("/","_")}.inf.lems${dat_suffix}.dat">
             <OutputColumn id="${g}" quantity="rampCellPop0[0]/test/${channel}/${g}/inf"/> 
         </OutputFile>
-        <OutputFile id="${g}_tau_of1" fileName="${channel}.${g}.tau.lems${dat_suffix}.dat">
+        <OutputFile id="${g}_tau_of1" fileName="${channel}.${g.replace("/","_")}.tau.lems${dat_suffix}.dat">
             <OutputColumn id="${g}" quantity="rampCellPop0[0]/test/${channel}/${g}/tau"/> 
         </OutputFile>
 #end## 
diff --git a/pyneuroml/analysis/NML2ChannelAnalysis.py b/pyneuroml/analysis/NML2ChannelAnalysis.py
index 7bb3c03c..66cd6f1f 100644
--- a/pyneuroml/analysis/NML2ChannelAnalysis.py
+++ b/pyneuroml/analysis/NML2ChannelAnalysis.py
@@ -469,9 +469,17 @@ def get_channel_gates(
         "gate_hh_rates",
         "gate_hh_tau_infs",
         "gate_hh_instantaneouses",
+        "gate_fractionals",
     ]:
         if hasattr(channel, gates):
-            channel_gates += [g.id for g in getattr(channel, gates)]
+            if gates == "gate_fractionals":
+                for g in getattr(channel, gates):
+                    for sg in g.sub_gates:
+                        channel_gates.append(str('%s/%s'%(g.id, sg.id)))
+            else:
+                for g in getattr(channel, gates):
+                    channel_gates.append(g.id)
+    #print('- Found gates: %s'%channel_gates)
     return channel_gates
 
 
diff --git a/pyneuroml/pynml.py b/pyneuroml/pynml.py
index af2fc3ee..560f52a0 100644
--- a/pyneuroml/pynml.py
+++ b/pyneuroml/pynml.py
@@ -897,7 +897,7 @@ def reload_standard_dat_file(file_name: str) -> typing.Tuple[dict, list]:
     """
     with open(file_name) as dat_file:
         data = {}  # type: dict
-        indeces = []  # type: list
+        indices = []  # type: list
         for line in dat_file:
             words = line.split()
 
@@ -905,13 +905,13 @@ def reload_standard_dat_file(file_name: str) -> typing.Tuple[dict, list]:
                 data["t"] = []
                 for i in range(len(words) - 1):
                     data[i] = []
-                    indeces.append(i)
+                    indices.append(i)
             data["t"].append(float(words[0]))
             for i in range(len(words) - 1):
                 data[i].append(float(words[i + 1]))
 
-        logger.info("Loaded data from %s; columns: %s" % (file_name, indeces))
-    return data, indeces
+        logger.info("Loaded data from %s; columns: %s" % (file_name, indices))
+    return data, indices
 
 
 def evaluate_component(comp_type, req_variables={}, parameter_values={}):
diff --git a/pyneuroml/utils/__init__.py b/pyneuroml/utils/__init__.py
index 85dbb784..3ced389b 100644
--- a/pyneuroml/utils/__init__.py
+++ b/pyneuroml/utils/__init__.py
@@ -327,6 +327,11 @@ def get_state_color(s: str) -> str:
     if s.startswith("u"):
         col = "#880088"
 
+    if '/' in s: # e.g. for sub gates
+        from pyneuroml.utils.plot import get_next_hex_color
+        col = get_next_hex_color()
+
+
     return col
 
 
diff --git a/pyneuroml/xppaut/__init__.py b/pyneuroml/xppaut/__init__.py
index 6c9f4822..e0b58f0e 100644
--- a/pyneuroml/xppaut/__init__.py
+++ b/pyneuroml/xppaut/__init__.py
@@ -6,6 +6,7 @@
 import typing
 import argparse
 import logging
+import os
 
 XPP_TIME = 'xpp_time'
 LEMS_TIME = 't'
@@ -79,6 +80,7 @@ def parse_script(file_path):
             elif line.startswith(('number', 'p', 'par', 'param')):
                 params = line.replace('number ', '').replace('par ', '').replace('p ', '').replace('param ', '')  # Skip the first word ('number', 'p', 'param' or 'par')
                 for pp1 in params.split(','):
+                    pp1 = pp1.replace(' = ','=').replace('= ','=').replace(' =','=') # better solution required...
                     for pp2 in pp1.split(' '):
                         if len(pp2)>0:
                             key, value = pp2.split('=')
@@ -189,7 +191,7 @@ def to_xpp(data, new_xpp_filename):
     xpp_ode+='\n# Parameters\n'
     for k,v in data["parameters"].items():
         if not k in INBUILT.keys():
-            xpp_ode += f'{k} = {v}\n'
+            xpp_ode += f'par {k} = {v}\n'
 
     xpp_ode+='\n# Functions\n'
     for k,v in data["functions"].items():
@@ -592,6 +594,79 @@ def main(args=None):
 
     cli(a=args)
 
+def run_xpp_file(filename, plot, show_plot_already=True, plot_separately={}):
+    import subprocess as sp
+    cmds = ['%s/xppaut'%os.environ["XPP_HOME"],filename, '-silent']
+    cwd = os.getcwd()
+    try:
+        ret_string = sp.check_output(
+            cmds, cwd=cwd, shell=False, stderr=sp.STDOUT
+        )
+        logger.info(
+            "Commands: %s completed successfully" % (cmds)
+        )
+        if isinstance(ret_string, bytes):
+            ret_string = ret_string.decode("utf-8")  # For Python 3...
+
+    except sp.CalledProcessError as err:
+        logger.error(
+            "CalledProcessError running commands: %s in %s (return code: %s), output:\n%s"
+            % (cmds, cwd, err.returncode, err.output),
+        )
+        raise err
+    except Exception as err:
+        logger.info(
+            "Error running commands: %s in (%s)!" % (cmds, cwd)
+        )
+        logger.info("Error: %s" % (err))
+        raise err
+    
+    if plot:
+
+        parsed_data = parse_script(filename)
+
+        from pyneuroml.pynml import reload_standard_dat_file
+        result_file = 'output.dat'
+        data, indices = reload_standard_dat_file(result_file)
+        logger.info('Loading data: %s'%(data.keys()))
+        default_figure = "Data loaded from: %s after running %s"%(result_file, filename)  # Title
+        ts = {default_figure:[]}
+        xs = {default_figure:[]}
+        labels = {default_figure:[]}
+        for new_fig in plot_separately:
+            ts[new_fig] = []
+            xs[new_fig] = []
+            labels[new_fig] = []
+
+        tds = list(parsed_data['time_derivatives'].keys())
+        cdvs = ['%s??'%c for c in parsed_data['conditional_derived_variables'].keys()]
+        outputs = tds+cdvs
+        logger.info('Loading data: %s, assuming these represent %s (%i values)'%(data.keys(),outputs, len(outputs)))
+        for i in indices:
+            label = outputs[i] if i < len(outputs) else '???'
+            
+            fig = default_figure
+            for new_fig in plot_separately:
+                if label in plot_separately[new_fig]:
+                    fig = new_fig
+            ts[fig].append(data['t'])
+            xs[fig].append(data[i])
+            labels[fig].append(label)
+
+        from pyneuroml.plot.Plot import generate_plot
+        axes = {}
+        for fig_title in ts:
+
+            ax = generate_plot(ts[fig_title], xs[fig_title],
+                title = fig_title,
+                labels=labels[fig_title],
+                xaxis="Time (?)",  # x axis legend
+                yaxis="??",  # y axis legend
+                show_plot_already=show_plot_already,  # Show or wait for plt.show()?
+            )  
+            axes[fig_title] = ax
+        return axes
+
 
 def cli(a: typing.Optional[typing.Any] = None, **kwargs: str):
     """Main cli caller method"""
@@ -608,19 +683,29 @@ def cli(a: typing.Optional[typing.Any] = None, **kwargs: str):
     lems_model_id = file_path.replace('.ode','').split('/')[-1]
     lems_model_file = file_path.replace('.ode','.model.xml')
 
+    if a.run and not a.lems and not a.xpp and not a.brian2:
+            
+        logger.info("Running %s with XPP (plotting: %s)..." % (a.ode_filename, a.plot))
+
+        run_xpp_file(a.ode_filename, a.plot)
+
+
     if a.lems:
         lems_filename = to_lems(parsed_data, lems_model_id, lems_model_file)
         logger.info("Generated LEMS file: %s" % lems_filename)
         if a.run:
             logger.info("Running %s with jNeuroML (plotting: %s)..." % (lems_filename, a.plot))
+
             from pyneuroml.runners import run_lems_with_jneuroml
             run_lems_with_jneuroml(lems_filename, nogui=True, plot=a.plot, verbose=True, load_saved_data=True)
 
 
     if a.xpp:
-        to_xpp(parsed_data, file_path.replace('.ode','_2.ode'))
-        if a.run or a.plot:
-            raise NotImplementedError("Running XPP not yet implemented!")
+        new_xppfile = file_path.replace('.ode','_2.ode')
+        to_xpp(parsed_data, new_xppfile)
+
+        if a.run:
+            run_xpp_file(new_xppfile, a.plot)
 
     if a.brian2:
         to_brian2(parsed_data, file_path.replace('.ode','_brian2.py'))