LEMS Export

Step-by-step guide to exporting TVBO models as LEMS XML

Introduction

This tutorial walks through the complete LEMS export workflow for three representative model types:

  1. A mean-field neural mass model (FitzHugh-Nagumo)
  2. A spiking integrate-and-fire model with spike/reset events
  3. A full simulation experiment with network, coupling, and output

1. FitzHugh-Nagumo — Basic Export

FitzHugh-Nagumo is already in the TVBO database as a Dynamics that maps cleanly to a LEMS ComponentType:

from tvbo import Dynamics
from tvbo.adapters.neuroml import NeuroMLAdapter

fhn = Dynamics.from_db("JansenRit")
adapter = NeuroMLAdapter(fhn)
xml = adapter.render_code()
print(xml)

<Lems>

  <!-- Tell jLEMS/jNeuroML which component is the simulation entry point. -->
  <Target component="sim_JansenRit"/>

  <!-- ════════════════════════════════════════════════════════════════
       Dimensions & Units (inline — no external includes needed)
       ════════════════════════════════════════════════════════════════ -->

  <!-- Dimensions -->
  <Dimension name="none"/>
  <Dimension name="time" t="1"/>
  <Dimension name="voltage" m="1" l="2" t="-3" i="-1"/>
  <Dimension name="per_time" t="-1"/>
  <Dimension name="conductance" m="-1" l="-2" t="3" i="2"/>
  <Dimension name="capacitance" m="-1" l="-2" t="4" i="2"/>
  <Dimension name="current" i="1"/>
  <Dimension name="resistance" m="1" l="2" t="-3" i="-2"/>
  <Dimension name="concentration" l="-3" n="1"/>
  <Dimension name="substance" n="1"/>
  <Dimension name="charge" t="1" i="1"/>
  <Dimension name="temperature" k="1"/>

  <!-- Units -->
  <Unit symbol="s" dimension="time" power="0"/>
  <Unit symbol="ms" dimension="time" power="-3"/>
  <Unit symbol="us" dimension="time" power="-6"/>
  <Unit symbol="V" dimension="voltage" power="0"/>
  <Unit symbol="mV" dimension="voltage" power="-3"/>
  <Unit symbol="A" dimension="current" power="0"/>
  <Unit symbol="mA" dimension="current" power="-3"/>
  <Unit symbol="nA" dimension="current" power="-9"/>
  <Unit symbol="pA" dimension="current" power="-12"/>
  <Unit symbol="S" dimension="conductance" power="0"/>
  <Unit symbol="mS" dimension="conductance" power="-3"/>
  <Unit symbol="nS" dimension="conductance" power="-9"/>
  <Unit symbol="F" dimension="capacitance" power="0"/>
  <Unit symbol="uF" dimension="capacitance" power="-6"/>
  <Unit symbol="nF" dimension="capacitance" power="-9"/>
  <Unit symbol="ohm" dimension="resistance" power="0"/>
  <Unit symbol="Mohm" dimension="resistance" power="6"/>
  <Unit symbol="per_s" dimension="per_time" power="0"/>
  <Unit symbol="per_ms" dimension="per_time" power="3"/>
  <Unit symbol="degC" dimension="temperature" offset="273.15"/>
  <Unit symbol="K" dimension="temperature" power="0"/>


  <!-- ════════════════════════════════════════════════════════════════
       LEMS Simulation infrastructure types
       (normally provided by Simulation.xml — defined inline to avoid
       the jNeuroML double-read bug)
       ════════════════════════════════════════════════════════════════ -->
  <ComponentType name="Simulation">
    <Parameter name="length" dimension="time"/>
    <Parameter name="step" dimension="time"/>
    <Children name="outputs" type="OutputFile"/>
    <ComponentReference name="target" type="Component"/>
    <Dynamics>
      <StateVariable name="t" dimension="time"/>
    </Dynamics>
    <Simulation>
      <Run component="target" variable="t" increment="step" total="length"/>
    </Simulation>
  </ComponentType>

  <ComponentType name="OutputFile">
    <Children name="outputColumn" type="OutputColumn"/>
    <Text name="fileName"/>
    <Text name="path"/>
    <Simulation>
      <DataWriter path="path" fileName="fileName"/>
    </Simulation>
  </ComponentType>

  <ComponentType name="OutputColumn">
    <Path name="quantity"/>
    <Simulation>
      <Record quantity="quantity"/>
    </Simulation>
  </ComponentType>

  <!-- ════════════════════════════════════════════════════════════════
       Dynamics ComponentType & Component instances
       ════════════════════════════════════════════════════════════════ -->

  <!-- ════════════════════════════════════════════════════════════════
       ComponentType: JansenRit
       Generated from TVBO Dynamics: JansenRit
       ════════════════════════════════════════════════════════════════ -->
  <ComponentType name="JansenRit">

    <!-- Parameters -->
    <Parameter name="A" dimension="none"/>
    <Parameter name="B" dimension="none"/>
    <Parameter name="J" dimension="none"/>
    <Parameter name="a_1" dimension="none"/>
    <Parameter name="a_2" dimension="none"/>
    <Parameter name="a_3" dimension="none"/>
    <Parameter name="a_4" dimension="none"/>
    <Parameter name="a" dimension="none"/>
    <Parameter name="b" dimension="none"/>
    <Parameter name="mu" dimension="none"/>
    <Parameter name="nu_max" dimension="none"/>
    <Parameter name="r" dimension="none"/>
    <Parameter name="v0" dimension="none"/>
    <!-- Coupling inputs -->
    <Parameter name="c_pop0" dimension="none"/>
    <Parameter name="local_coupling" dimension="none"/>
    <!-- Initial condition parameters -->
    <Parameter name="y0_0" dimension="none"/>
    <Parameter name="y1_0" dimension="none"/>
    <Parameter name="y2_0" dimension="none"/>
    <Parameter name="y3_0" dimension="none"/>
    <Parameter name="y4_0" dimension="none"/>
    <Parameter name="y5_0" dimension="none"/>

    <!-- Time constant for derivatives.
         All parameters use dimension="none" to avoid jLEMS dimensional
         analysis failures (neural-mass models mix units freely).
         / SEC is always applied so the TimeDerivative has per_time
         dimension.  When the model already has time units in parameters,
         SEC = 1s (numerically 1.0 — identity); otherwise SEC = 1<time_scale>
         provides the actual numeric conversion.
         needs_sec=True  time_scale=ms -->
    <Constant name="SEC" dimension="time" value="1ms"/>

    <!-- Exposures (one per state variable) -->
    <Exposure name="y0" dimension="none"/>
    <Exposure name="y1" dimension="none"/>
    <Exposure name="y2" dimension="none"/>
    <Exposure name="y3" dimension="none"/>
    <Exposure name="y4" dimension="none"/>
    <Exposure name="y5" dimension="none"/>

    <Dynamics>

      <!-- State variables -->
      <StateVariable name="y0" dimension="none" exposure="y0"/>
      <StateVariable name="y1" dimension="none" exposure="y1"/>
      <StateVariable name="y2" dimension="none" exposure="y2"/>
      <StateVariable name="y3" dimension="none" exposure="y3"/>
      <StateVariable name="y4" dimension="none" exposure="y4"/>
      <StateVariable name="y5" dimension="none" exposure="y5"/>

      <!-- Derived variables (simple and conditional/piecewise) -->
      <DerivedVariable name="sigma_y0_1" dimension="none" value="2.0*nu_max/(1.0 + exp(r*(v0 - J*a_1*y0)))"/>
      <DerivedVariable name="sigma_y0_3" dimension="none" value="2.0*nu_max/(1.0 + exp(r*(v0 - J*a_3*y0)))"/>
      <DerivedVariable name="sigma_y1_y2" dimension="none" value="2.0*nu_max/(1.0 + exp(r*(v0 + y2 - y1)))"/>

      <!-- ── Flat dynamics (no spike events) ── -->

      <!-- Time derivatives -->
      <TimeDerivative variable="y0" value="(y3) / SEC"/>
      <TimeDerivative variable="y1" value="(y4) / SEC"/>
      <TimeDerivative variable="y2" value="(y5) / SEC"/>
      <TimeDerivative variable="y3" value="(-y0*a^2 - 2.0*a*y3 + A*a*sigma_y1_y2) / SEC"/>
      <TimeDerivative variable="y4" value="(-y1*a^2 - 2.0*a*y4 + A*a*(c_pop0 + mu + local_coupling*(y1 - y2) + J*a_2*sigma_y0_1)) / SEC"/>
      <TimeDerivative variable="y5" value="(-y2*b^2 - 2.0*b*y5 + B*J*a_4*b*sigma_y0_3) / SEC"/>

      <!-- Initial conditions -->
      <OnStart>
        <StateAssignment variable="y0" value="y0_0"/>
        <StateAssignment variable="y1" value="y1_0"/>
        <StateAssignment variable="y2" value="y2_0"/>
        <StateAssignment variable="y3" value="y3_0"/>
        <StateAssignment variable="y4" value="y4_0"/>
        <StateAssignment variable="y5" value="y5_0"/>
      </OnStart>

      <!-- Events (non-spike) -->

    </Dynamics>

  </ComponentType>


  <!-- ════════════════════════════════════════════════════════════════
       Component instances (default parameter values)
       ════════════════════════════════════════════════════════════════ -->
  <Component id="JansenRit_inst" type="JansenRit" A="3.25" B="22.0" J="135.0" a_1="1.0" a_2="0.8" a_3="0.25" a_4="0.25" a="0.1" b="0.05" mu="0.22" nu_max="0.0025" r="0.56" v0="5.52" c_pop0="0" local_coupling="0" y0_0="0.1" y1_0="0.1" y2_0="0.1" y3_0="0.1" y4_0="0.1" y5_0="0.1"/>



  <!-- ════════════════════════════════════════════════════════════════
       Simulation — target the component instance directly.
       Output: one file with all state variables.
       ════════════════════════════════════════════════════════════════ -->
  <Simulation id="sim_JansenRit" length="1000.0ms" step="0.01220703125ms" target="JansenRit_inst">

    <OutputFile id="of1" fileName="results/JansenRit.dat">
      <OutputColumn id="y0" quantity="y0"/>
      <OutputColumn id="y1" quantity="y1"/>
      <OutputColumn id="y2" quantity="y2"/>
      <OutputColumn id="y3" quantity="y3"/>
      <OutputColumn id="y4" quantity="y4"/>
      <OutputColumn id="y5" quantity="y5"/>
    </OutputFile>

  </Simulation>

</Lems>

The rendered XML includes:

  • <ComponentType name="FitzHughNagumo"> with all model parameters as <Parameter>
  • <StateVariable name="V" ...> and <StateVariable name="W" ...> with <Exposure>
  • <TimeDerivative variable="V" value="(V - V^3/3 - W + I) / SEC"/> — note the / SEC to account for the TVBO convention of millisecond-based integration with a 1s constant
  • <OnStart> with <StateAssignment> from initial conditions
  • A <Component id="FitzHughNagumo" ...> instance with default parameter values
  • A <Simulation> block with output files per state variable

2. Integrate-and-Fire — Spike/Reset Events

Define an IaF model in TVBO YAML:

# database/dynamics/spiking/IaF.yaml
name: IaF
iri: "GO:0019228"  # neuronal action potential

parameters:
  tau_m:
    value: 20.0
    unit: ms
  V_th:
    value: -55.0
    unit: mV
  V_reset:
    value: -70.0
    unit: mV
  I_ext:
    value: 15.0

state_variables:
  V:
    equation:
      rhs: "(V_rest - V + I_ext) / tau_m"
    unit: mV
    initial_value: -70.0

events:
  spike:
    condition:
      rhs: "V >= V_th"
    affect:
      rhs: "V = V_reset"

The adapter maps events.spike.condition.rhs = "V >= V_th" to:

<OnCondition test="V .geq. V_th">
  <StateAssignment variable="V" value="V_reset"/>
  <EventOut port="spike"/>
</OnCondition>

The relational operators (>=.geq.) are handled by LEMSPrinter inside sympy_to_lems().


3. Conditional Derived Variables

Gate variables in Hodgkin-Huxley-style models use piecewise functions with a singularity at specific voltages. TVBO stores these as SymPy Piecewise expressions, which the template automatically maps to LEMS <ConditionalDerivedVariable>:

derived_variables:
  alpha_m:
    equation:
      rhs: "Piecewise((1.0, Eq(V, -40)), (0.1*(V+40)/(1-exp(-(V+40)/10)), True))"
    unit: per_ms

Exported as:

<ConditionalDerivedVariable name="alpha_m" dimension="per_time">
  <Case condition="V .eq. -40" value="1.0"/>
  <Case value="0.1*(V + 40)/(1 - exp(-(V + 40)/10))"/>
</ConditionalDerivedVariable>

The final <Case> with no condition is the default fallback — equivalent to True in SymPy.


4. Coupling

For a whole-brain experiment, TVBO’s Coupling metadata (pre-expression, post-expression, global coupling constant) is exported as a dedicated LEMS ComponentType:

<ComponentType name="Coupling">
  <Parameter name="global_coupling" dimension="none"/>
  <Dynamics>
    <DerivedVariable name="pre"  dimension="none" value="V_j"/>
    <DerivedVariable name="gx"   dimension="none" value="global_coupling * pre"/>
    <DerivedVariable name="post" dimension="none" value="a * gx"/>
    <DerivedParameter name="c_pop0" value="post"/>
  </Dynamics>
</ComponentType>

The c_pop0 DerivedParameter is the coupling term that appears in the network node’s state variable equations.


5. Full Experiment Export

from tvbo import SimulationExperiment
from tvbo.adapters.neuroml import NeuroMLAdapter

exp = SimulationExperiment.from_string("""
label: JansenRit LEMS Export
dynamics:
  name: JansenRit
""")
adapter = NeuroMLAdapter(exp)

# Render monolithic LEMS file
xml = adapter.render_code()

# Validate with PyLEMS
try:
    adapter.validate(xml)
    print("Valid LEMS")
except Exception as e:
    print(f"Validation error: {e}")

# Save to file
with open("lems_sim.xml", "w") as f:
    f.write(xml)
Valid LEMS

6. Split-File Export

By NeuroML convention, a fully portable model is distributed as three separate LEMS files that references each other via <Include>. The TVBO adapter generates all three with a single call:

paths = adapter.export("./output", split=True, validate=False)
print(paths)
# {
#   'dynamics':   './output/ses-..._dynamics.xml',
#   'network':    './output/ses-..._network.xml',
#   'simulation': './output/ses-..._simulation.xml',
# }
{'dynamics': 'output/ses-2_desc-None_dynamics.xml', 'network': 'output/ses-2_desc-None_network.xml', 'simulation': 'output/ses-2_desc-None_simulation.xml'}

The chain works as follows:

File Contains Includes
*_dynamics.xml Dimensions, Units, ComponentType, Coupling CT, default instances nothing
*_network.xml <Component id="net" ...> + population *_dynamics.xml
*_simulation.xml <Simulation> + Output files *_network.xml

You can render each piece independently:

# Just the ComponentType definitions
xml_dyn = adapter.render_dynamics()

# Network file that references an external dynamics file
xml_net = adapter.render_network(dynamics_file="my_dynamics.xml")

# Simulation file that references an external network file
xml_sim = adapter.render_simulation(network_file="my_network.xml")

7. Unified render() API

SimulationExperiment.render() is the unified entry point for all serialisation formats — both executable code and metadata.

# exp and adapter are already defined in section 5 above

# Render as LEMS XML (same as NeuroMLAdapter(exp).render_code())
xml = exp.render("lems")       # also accepts 'neuroml', 'nml'

# Render as YAML specification
yaml_str = exp.render("yaml")

# Render as JAX code
jax_code = exp.render("jax")

For split-file export or other adapter-specific options, use NeuroMLAdapter directly.


8. Running via jNeuroML

Once a valid LEMS file is produced, simulate it with jNeuroML:

# Via the adapter
adapter.run()

# Via SimulationExperiment.run()
exp.run("neuroml")
JansenRit LEMS Export
└── integration
        data: (81921, 6)

Or from the command line:

jnml ses-..._simulation.xml

Results are written to results/<model>_<sv>.dat (one file per state variable per run), matching the <OutputFile> definitions in the <Simulation> block.


9. User-Defined Functions

LEMS has no mechanism for user-defined functions. Models that define helper functions in YAML (e.g. a sigmoid Sigm(x)) have those calls automatically inlined by the adapter before rendering:

# In TVBO YAML — user function
functions:
  Sigm:
    arguments: [x]
    equation:
      rhs: "2*e0 / (1 + exp(r*(v0 - x)))"

state_variables:
  y0:
    equation:
      rhs: "A*a * Sigm(y1 - y2) - 2*a*y3 - a**2*y0"

The adapter substitutes Sigm(y1 - y2) with its body, producing valid LEMS:

<TimeDerivative variable="y0"
  value="(A*a * (2*e0/(1 + exp(r*(v0 - (y1 - y2))))) - 2*a*y3 - a^2*y0) / SEC"/>

This inlining is done symbolically via SymPy’s replace() + xreplace(), so nested or composed function calls work correctly.


10. Understanding the LEMS Structure

A complete TVBO → LEMS export produces this XML skeleton:

<Lems>
  <!-- Standard NeuroML dimensions (voltage, time, current, ...) -->
  <Dimension name="voltage" m="1" l="2" t="-3" i="-1"/>
  ...

  <!-- Units (mV, ms, nA, ...) -->
  <Unit name="milliVolt" symbol="mV" dimension="voltage" power="-3"/>
  ...

  <!-- The dynamics model as a ComponentType -->
  <ComponentType name="MyModel">
    <Parameter name="param1" dimension="none"/>
    ...
    <Constant name="SEC" dimension="time" value="1s"/>
    <Exposure name="V" dimension="voltage"/>

    <Dynamics>
      <StateVariable name="V" dimension="voltage" exposure="V"/>

      <!-- Simple derived variable -->
      <DerivedVariable name="f_V" dimension="none" value="1/(1 + exp(-V))"/>

      <!-- Piecewise derived variable -->
      <ConditionalDerivedVariable name="alpha_n" dimension="per_time">
        <Case condition="V .eq. -55" value="0.1"/>
        <Case value="0.01*(V + 55)/(1 - exp(-(V+55)/10))"/>
      </ConditionalDerivedVariable>

      <!-- Time derivative (scaled to ms time base) -->
      <TimeDerivative variable="V" value="(dV_rhs) / SEC"/>

      <!-- Initial conditions -->
      <OnStart>
        <StateAssignment variable="V" value="V_0"/>
      </OnStart>

      <!-- Spike/reset events -->
      <OnCondition test="V .geq. V_th">
        <StateAssignment variable="V" value="V_reset"/>
      </OnCondition>
    </Dynamics>
  </ComponentType>

  <!-- Coupling function as a separate ComponentType -->
  <ComponentType name="Coupling">
    ...
  </ComponentType>

  <!-- Component instance with default parameter values -->
  <Component id="MyModel" type="MyModel" param1="1.0" V_0="-65.0"/>

  <!-- Network: population of N nodes -->
  <Component id="net" type="network">
    <Component id="pop0" type="population" component="MyModel" size="80"/>
  </Component>

  <!-- Simulation: length, step, output files -->
  <Simulation id="sim1" length="1000ms" step="0.1ms" target="net">
    <OutputFile id="of_V" fileName="results/MyModel_V.dat">
      <OutputColumn id="V_0" quantity="pop0[0]/V"/>
      ...
    </OutputFile>
  </Simulation>
</Lems>

Troubleshooting

DeprecationWarning: Dynamics.to_lems()
This method is deprecated. Use NeuroMLAdapter(model).render_code() instead.
LEMSPrinter output has exp() instead of e^x
That’s correct — LEMS uses the function name exp(x), not e^x.
Validation fails with “Unknown element: network”
The TVBO template generates self-contained LEMS, not standard NeuroML XML. It does not use <neuroml> as the root element — <Lems> is the root. Use jnml directly, not pynml_run_neuroml2_lems.
Missing units dimension
If a parameter has unit: "mA/cm2" that isn’t in the lookup table, the dimension falls back to "none". Extend UNITS in tvbo/adapters/neuroml.py or map it in your YAML with a known unit string.