diff --git a/bmtk/simulator/pointnet/pointnetwork.py b/bmtk/simulator/pointnet/pointnetwork.py
index 89877b351..e3852caed 100644
--- a/bmtk/simulator/pointnet/pointnetwork.py
+++ b/bmtk/simulator/pointnet/pointnetwork.py
@@ -76,7 +76,6 @@ def __init__(self, **properties):
         self._params_cache = {}
 
         self._virtual_ids_map = {}
-
         self._batch_nodes = True
 
         # self._nest_id_map = {}
@@ -90,6 +89,8 @@ def __init__(self, **properties):
         self._gid_map = GidPool()
         self._virtual_gids = GidPool()
 
+        self._nestml_models = []
+
     @property
     def py_function_caches(self):
         return pyfunction_cache
@@ -168,6 +169,9 @@ def build_recurrent_edges(self, force_resolution=False):
             for edge in edge_pop.get_edges():
                 nest_srcs = self.gid_map.get_nestids(edge_pop.source_nodes, edge.source_node_ids)
                 nest_trgs = self.gid_map.get_nestids(edge_pop.target_nodes, edge.target_node_ids)
+                if np.isscalar(edge.nest_params['weight']):
+                    edge.nest_params['weight'] = np.full(shape=len(nest_srcs),
+                                                         fill_value=edge.nest_params['weight'])
                 self._nest_connect(nest_srcs, nest_trgs, conn_spec='one_to_one', syn_spec=edge.nest_params)
 
     def find_edges(self, source_nodes=None, target_nodes=None):
@@ -219,6 +223,9 @@ def add_spike_trains(self, spike_trains, node_set, sg_params={'precise_times': T
                 for edge in edge_pop.get_edges():
                     nest_trgs = self.gid_map.get_nestids(edge_pop.target_nodes, edge.target_node_ids)
                     nest_srcs = virt_gid_map.get_nestids(edge_pop.source_nodes, edge.source_node_ids)
+                    if np.isscalar(edge.nest_params['weight']):
+                        edge.nest_params['weight'] = np.full(shape=len(nest_srcs),
+                                                             fill_value=edge.nest_params['weight'])
                     self._nest_connect(nest_srcs, nest_trgs, conn_spec='one_to_one', syn_spec=edge.nest_params)
 
     def _nest_connect(self, nest_srcs, nest_trgs, conn_spec='one_to_one', syn_spec=None):
@@ -230,7 +237,7 @@ def _nest_connect(self, nest_srcs, nest_trgs, conn_spec='one_to_one', syn_spec=N
             # An occuring issue is when dt > delay, add some extra messaging in log to help users fix problem.
             res_kernel = nest.GetKernelStatus().get('resolution', 'NaN')
             delay_edges = syn_spec.get('delay', 'NaN')
-            msg = 'synaptic "delay" value in edges ({}) is not compatable with simulator resolution/"dt" ({})'.format(
+            msg = 'synaptic "delay" value in edges ({}) is not compatible with simulator resolution/"dt" ({})'.format(
                 delay_edges, res_kernel
             )
             self.io.log_error('{}{}'.format(bde.errorname, bde.errormessage))
@@ -241,3 +248,25 @@ def _nest_connect(self, nest_srcs, nest_trgs, conn_spec='one_to_one', syn_spec=N
             # Record exception to log file.
             self.io.log_error(str(e))
             raise
+
+    def add_nestml_models(self, modelname):
+        if modelname not in self._nestml_models:
+            self._nestml_models.append(modelname)
+
+    def initialize_nestml(self, rebuild_nestml=True):
+        from pynestml.frontend.pynestml_frontend import generate_nest_target
+        models_dir = os.path.commonpath([self.get_component('point_neuron_models_dir'), self.get_component('synaptic_models_dir')])
+        nestml_model_name_suffix = '_nestml'
+
+        if self._nestml_models:
+            if rebuild_nestml or not (os.path.exists(os.path.join(models_dir, 'nestml_target'))):
+                if not rebuild_nestml:
+                    self.io.log_info('Cannot find nestml target, must rebuild')
+                generate_nest_target(input_path=self._nestml_models,
+                                     logging_level='ERROR',
+                                     suffix=nestml_model_name_suffix,
+                                     target_path=os.path.join(models_dir, 'nestml_target'))
+                self.io.log_info('Generating NESTML models: ', self._nestml_models)
+                # a module by the same name can only be loaded once; do this at the very end of the function
+            nest.Install('nestmlmodule')
+
diff --git a/bmtk/simulator/pointnet/pointsimulator.py b/bmtk/simulator/pointnet/pointsimulator.py
index bd5f27ece..3c1454e9b 100644
--- a/bmtk/simulator/pointnet/pointsimulator.py
+++ b/bmtk/simulator/pointnet/pointsimulator.py
@@ -47,6 +47,7 @@ def __init__(self, graph, dt=0.001, overwrite=True, print_time=False):
 
         self._cells_built = False
         self._internal_connections_built = False
+        self._rebuild_nestml = True
 
         self._graph = graph
         self._external_cells = {}  # dict-of-dict of external pointnet cells with keys [network_name][cell_id]
@@ -228,6 +229,11 @@ def from_config(cls, configure, graph):
         if run_dict.get('allow_offgrid_spikes', False):
             network.set_spike_generator_params(allow_offgrid_spikes=True)
 
+        if 'rebuild_nestml' in config:
+            network._rebuild_nestml = config['rebuild_nestml']
+
+        graph.initialize_nestml(network._rebuild_nestml)
+
         # Create the output-directory, or delete existing files if it already exists
         graph.io.log_info('Setting up output directory')
         if not os.path.exists(config['output']['output_dir']):
@@ -276,7 +282,7 @@ def from_config(cls, configure, graph):
                 mod = mods.SpikesMod(**report.params)
 
             elif isinstance(report, reports.MembraneReport):
-                # For convience and for compliance with SONATA format. "membrane_report" and "multimeter_report is the
+                # For convenience and for compliance with SONATA format. "membrane_report" and "multimeter_report is the
                 # same in pointnet.
                 mod = mods.MultimeterMod(**report.params)
 
diff --git a/bmtk/simulator/pointnet/sonata_adaptors.py b/bmtk/simulator/pointnet/sonata_adaptors.py
index bb3296f95..342cb08e1 100644
--- a/bmtk/simulator/pointnet/sonata_adaptors.py
+++ b/bmtk/simulator/pointnet/sonata_adaptors.py
@@ -1,7 +1,9 @@
+import logging
 import numpy as np
 from collections import Counter
 import numbers
 import nest
+import os
 import types
 import pandas as pd
 
@@ -11,6 +13,8 @@
 from bmtk.simulator.pointnet.glif_utils import convert_aibs2nest
 from bmtk.simulator.pointnet.nest_utils import nest_version
 
+logger = logging.getLogger(__name__)
+
 NEST_SYNAPSE_MODEL_PROP = 'model' if nest_version[0] == 2 else 'synapse_model'
 
 
@@ -155,6 +159,7 @@ def get_batches(self, node_group):
     def preprocess_node_types(network, node_population):
         NodeAdaptor.preprocess_node_types(network, node_population)
         node_types_table = node_population.types_table
+
         if 'model_template' in node_types_table.columns and 'dynamics_params' in node_types_table.columns:
             node_type_ids = np.unique(node_population.type_ids)
             for nt_id in node_type_ids:
@@ -166,6 +171,20 @@ def preprocess_node_types(network, node_population):
                     node_type_attrs['model_template'] = model_template
                     node_type_attrs['dynamics_params'] = dynamics_params
 
+                if mtemplate.startswith('nestml:'):
+                    # NESTML model was requested
+
+                    model_name = mtemplate.split(':')[1]
+
+                    models_dir = network.get_component('point_neuron_models_dir')
+                    model_fn = os.path.join(models_dir, model_name + '.nestml')
+                    network.add_nestml_models(model_fn)
+
+                    # replace model_template name with suffixed version
+
+                    node_types_table[nt_id]['model_template'] = 'nest:' + model_name + '_nestml'
+                    node_types_table._df_cache = None   # clear the table internal cache
+
     @staticmethod
     def patch_adaptor(adaptor, node_group, network):
         node_adaptor = NodeAdaptor.patch_adaptor(adaptor, node_group, network)
@@ -256,14 +275,24 @@ def preprocess_edge_types(network, edge_population):
         # Fix for sonata/300_pointneurons
         EdgeAdaptor.preprocess_edge_types(network, edge_population)
         edge_types_table = edge_population.types_table
-        edge_type_ids = np.unique(edge_population.type_ids)
-
-        for et_id in edge_type_ids:
-            edge_type = edge_types_table[et_id]
-            if 'model_template' in edge_types_table.columns:
-                model_template = edge_type['model_template']
-                if model_template.startswith('nest'):
-                    edge_type['model_template'] = model_template[5:]
+
+        if 'model_template' in edge_types_table.columns and 'dynamics_params' in edge_types_table.columns:
+            edge_type_ids = np.unique(edge_population.type_ids)
+            for et_id in edge_type_ids:
+                edge_type_attrs = edge_types_table[et_id]
+                mtemplate = edge_type_attrs['model_template']
+                if mtemplate.startswith('nest:'):
+                    edge_type['model_template'] = mtemplate[5:]
+
+                if mtemplate.startswith('nestml:'):
+                    model_name = mtemplate.split(':')[1]
+
+                    models_dir = network.get_component('synaptic_models_dir')
+                    model_fn = os.path.join(models_dir, model_name + '.nestml')
+                    network.add_nestml_models(model_fn)
+
+                    edge_types_table[et_id]['model_template'] = model_name + '_nestml'
+                    edge_types_table._df_cache = None  # clear the table internal cache
 
     def get_batches(self, edge_group):
         src_ids = {}
diff --git a/examples/point_120cells_nestml/README.md b/examples/point_120cells_nestml/README.md
new file mode 100644
index 000000000..40de8f7f4
--- /dev/null
+++ b/examples/point_120cells_nestml/README.md
@@ -0,0 +1,69 @@
+# 120 point neuron network
+
+A small network of 120 point-neurons. Neuron and synapse models are specified as NESTML models to demonstrate how to incorporate NESTML models. Uses PointNet and will require NEST and NESTML to run.
+
+
+## Running:
+To run a simulation, install bmtk and run the following:
+```
+$ python run_pointnet.py config.simulation.json
+```
+If successful, will create a *output* directory containing log, spike trains and recorded cell variables.
+
+## The Network:
+The network files have already been built and stored as SONATA files in the *network/* directory. The bmtk Builder
+script used to create the network files is *build_network.py*. To adjust the parameters and/or topology of the network
+change this file and run:
+```
+$ python build_network.py
+```
+This will overwrite the existing files in the network directory. Note that there is some randomness in how the network
+is built, so expect (slightly) different simulation results everytime the network is rebuilt.
+
+## Simulation Parameters
+Parameters to run the simulation, including run-time, inputs, recorded variables, and networks are stored in config.json
+and can modified with a text editor.
+
+## Plotting results
+```
+$ python plot_output.py
+```
+
+## Perturbation simulations
+The file ```config_pertubration.json``` uses the same network (recurrent connections + feedforward inputs) as before. However
+is designed to also simulate the effects of optogenetic or current clamp inhibition and excitation on a subset of the cells. 
+To run this example:
+```bash
+$ python run_pointnet.py config_perturbations.json
+``` 
+
+The only difference between this simulation and the original is in the **inputs** section in the config. To simulate the 
+perturbations we use a current clamp input to make cells 20 through 39 overly excitatory:
+```json
+"exc_perturbation": {
+   "input_type": "current_clamp",
+   "module": "IClamp",
+   "node_set": {
+     "population": "cortex",
+     "node_ids": [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]
+   },
+   "amp": 230.0,
+   "delay": 1.0,
+   "duration": 3000.0
+},
+```
+
+For cells 40 through 49 we can use a negative current to depress the cell activity
+```json
+"inh_perturbation": {
+  "input_type": "current_clamp",
+  "module": "IClamp",
+  "node_set": {
+    "population": "cortex",
+    "node_ids": [40, 41, 42, 43, 44, 45, 46, 47, 48, 49]
+  },
+  "amp": -230.0,
+  "delay": 1.0,
+  "duration": 3000.0
+}
+```
diff --git a/examples/point_120cells_nestml/build_network.py b/examples/point_120cells_nestml/build_network.py
new file mode 100644
index 000000000..568d513bf
--- /dev/null
+++ b/examples/point_120cells_nestml/build_network.py
@@ -0,0 +1,116 @@
+import numpy as np
+
+from bmtk.builder import NetworkBuilder
+from bmtk.builder.auxi.node_params import positions_columinar
+
+
+lif_models = {
+    'LIF_exc': {
+        'N': 80,
+        'ei': 'e',
+        'pop_name': 'LIF_exc',
+        'model_type': 'point_neuron',
+        'model_template': 'nestml:iaf_psc_alpha',
+        'dynamics_params': 'iaf_psc_delta_exc.json',
+    },
+    'LIF_inh': {
+        'N': 40,
+        'ei': 'i',
+        'pop_name': 'LIF_inh',
+        'model_type': 'point_neuron',
+        'model_template': 'nestml:iaf_psc_delta',
+        'dynamics_params': 'iaf_psc_delta_inh.json',
+    }
+}
+
+input_network_model = {
+    'input_network': {
+        'N': 100,
+        'ei': 'e',
+        'pop_name': 'input_network',
+        'model_type': 'virtual'
+    }
+}
+
+
+def random_connections(source, target, p=0.1):
+    sid = source['node_id']  # Get source id
+    tid = target['node_id']  # Get target id
+
+    # Avoid self-connections.
+    if sid == tid:
+        return None
+
+    return np.random.binomial(1, p)  # nsyns
+
+
+def build_cortex_network():
+    cortex = NetworkBuilder('cortex')
+    for model in lif_models:
+        params = lif_models[model].copy()
+        positions = positions_columinar(N=lif_models[model]['N'], center=[0, 10.0, 0], max_radius=50.0, height=200.0)
+        cortex.add_nodes(
+            x=positions[:, 0],
+            y=positions[:, 1],
+            z=positions[:, 2],
+            **params
+        )
+
+    cortex.add_edges(
+        source={'ei': 'e'},
+        connection_rule=random_connections,
+        connection_params={'p': 0.1},
+        syn_weight=2.0,  # 2.0
+        delay=1.5,
+        dynamics_params='ExcToInh.json',
+        model_template='static_synapse'
+    )
+
+    cortex.add_edges(
+        source={'ei': 'i'},
+        connection_rule=random_connections,
+        connection_params={'p': 0.1},
+        syn_weight=-1.5,
+        delay=1.5,
+        dynamics_params='InhToExc.json',
+        model_template='static_synapse'
+    )
+
+    cortex.build()
+    cortex.save(output_dir='network')
+
+    return cortex
+
+
+def build_thalamus_network(cortex):
+    thalamus = NetworkBuilder('thalamus')
+    thalamus.add_nodes(**input_network_model['input_network'])
+
+    thalamus.add_edges(
+        target=cortex.nodes(ei='e'),
+        connection_rule=random_connections,
+        connection_params={'p': 0.1},
+        syn_weight=220.0, 
+        delay=1.5,
+        dynamics_params='ExcToExc.json',
+        model_template='nestml:static_synapse'
+    )
+
+    thalamus.add_edges(
+        target=cortex.nodes(ei='i'),
+        connection_rule=random_connections,
+        connection_params={'p': 0.1},
+        syn_weight=5.0,
+        delay=1.5,
+        dynamics_params='ExcToExc.json',
+        model_template='nestml:static_synapse'
+    )
+    thalamus.build()
+    thalamus.save(output_dir='network')
+
+    return thalamus
+
+
+if __name__ == '__main__':
+    cortex_net = build_cortex_network()
+    thalamus_net = build_thalamus_network(cortex_net)
diff --git a/examples/point_120cells_nestml/config.circuit.json b/examples/point_120cells_nestml/config.circuit.json
new file mode 100644
index 000000000..fa6f70e08
--- /dev/null
+++ b/examples/point_120cells_nestml/config.circuit.json
@@ -0,0 +1,36 @@
+{
+  "manifest": {
+	"$BASE_DIR": ".",
+    "$NETWORK_DIR": "$BASE_DIR/../point_120cells/network",
+    "$LOCAL_NETWORK_DIR": "$BASE_DIR/network",
+    "$MODELS_DIR": "$BASE_DIR/../point_components"
+  },
+
+  "components": {
+    "point_neuron_models_dir": "$MODELS_DIR/cell_models",
+	"synaptic_models_dir": "$MODELS_DIR/synaptic_models"
+  },
+
+  "networks": {
+    "nodes": [
+      {
+        "nodes_file": "$NETWORK_DIR/cortex_nodes.h5",
+        "node_types_file": "$LOCAL_NETWORK_DIR/cortex_node_types.csv"
+      },
+      {
+        "nodes_file": "$NETWORK_DIR/thalamus_nodes.h5",
+        "node_types_file": "$NETWORK_DIR/thalamus_node_types.csv"
+      }
+    ],
+    "edges": [
+      {
+        "edges_file": "$NETWORK_DIR/cortex_cortex_edges.h5",
+        "edge_types_file": "$NETWORK_DIR/cortex_cortex_edge_types.csv"
+      },
+      {
+        "edges_file": "$NETWORK_DIR/thalamus_cortex_edges.h5",
+        "edge_types_file": "$NETWORK_DIR/thalamus_cortex_edge_types.csv"
+      }
+    ]
+  }
+}
diff --git a/examples/point_120cells_nestml/config.simulation.json b/examples/point_120cells_nestml/config.simulation.json
new file mode 100644
index 000000000..9fb2c3b8b
--- /dev/null
+++ b/examples/point_120cells_nestml/config.simulation.json
@@ -0,0 +1,50 @@
+{
+  "manifest": {
+	"$BASE_DIR": "${configdir}",
+    "$NETWORK_DIR": "$BASE_DIR/network",
+    "$MODELS_DIR": "$BASE_DIR/../point_components",
+    "$OUTPUT_DIR": "$BASE_DIR/output",
+    "$INPUT_DIR": "$BASE_DIR/inputs"
+  },
+
+  "run": {
+    "tstop": 3000.0,
+    "dt": 0.001,
+    "block_run": false,
+    "block_size": 1000.0
+  },
+
+  "inputs": {
+    "thalamus_spikes": {
+      "input_type": "spikes",
+      "module": "sonata",
+      "input_file": "$INPUT_DIR/thalamus_spikes.h5",
+      "node_set": "thalamus"
+    }
+  },
+
+  "reports": {
+    "membrane_potential": {
+      "cells": {
+        "population": "cortex",
+        "node_id": [0, 20, 60, 80, 100]
+      },
+      "variable_name": "V_m",
+      "module": "multimeter_report",
+      "sections": "soma"
+    }
+  },
+
+  "output": {
+	"log_file": "log.txt",
+	"spikes_file": "spikes.h5",
+    "spikes_file_csv": "spikes.csv",
+	"output_dir": "$OUTPUT_DIR",
+    "overwrite_output_dir": true,
+    "quiet_simulator": true
+  },
+
+  "target_simulator":"NEST",
+  "rebuild_nestml": true,
+  "network": "config.circuit.json"
+}
diff --git a/examples/point_120cells_nestml/config.simulation_perturbations.json b/examples/point_120cells_nestml/config.simulation_perturbations.json
new file mode 100644
index 000000000..9fc38daa9
--- /dev/null
+++ b/examples/point_120cells_nestml/config.simulation_perturbations.json
@@ -0,0 +1,69 @@
+{
+  "manifest": {
+	"$BASE_DIR": ".",
+    "$OUTPUT_DIR": "$BASE_DIR/output",
+    "$INPUT_DIR": "$BASE_DIR/inputs"
+  },
+
+  "target_simulator":"NEST",
+  "rebuild_nestml":true,
+  "run": {
+    "tstop": 3000.0,
+    "dt": 0.001
+  },
+
+  "inputs": {
+    "LGN_spikes": {
+      "input_type": "spikes",
+      "module": "csv",
+      "input_file": "thalamus_spikes.csv",
+      "node_set": "thalamus"
+    },
+
+    "exc_perturbation": {
+      "input_type": "current_clamp",
+      "module": "IClamp",
+      "node_set": {
+        "population": "cortex",
+        "node_ids": [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]
+      },
+      "amp": 230.0,
+      "delay": 1.0,
+      "duration": 3000.0
+    },
+
+    "inh_perturbation": {
+      "input_type": "current_clamp",
+      "module": "IClamp",
+      "node_set": {
+        "population": "cortex",
+        "node_ids": [40, 41, 42, 43, 44, 45, 46, 47, 48, 49]
+      },
+      "amp": -230.0,
+      "delay": 1.0,
+      "duration": 3000.0
+    }
+
+  },
+
+  "reports": {
+    "membrane_potential": {
+      "cells": {"population": "cortex", "node_id": [0, 20, 60, 80, 100]},
+      "variable_name": "V_m",
+      "module": "multimeter_report",
+      "sections": "soma",
+      "enabled": true
+    }
+  },
+
+  "output": {
+	"log_file": "log.txt",
+	"spikes_file": "spikes.h5",
+    "spikes_file_csv": "spikes.csv",
+	"output_dir": "$OUTPUT_DIR",
+    "overwrite_output_dir": true,
+    "quiet_simulator": true
+  },
+
+  "network": "config.circuit.json"
+}
diff --git a/examples/point_120cells_nestml/create_inputs.py b/examples/point_120cells_nestml/create_inputs.py
new file mode 100644
index 000000000..058433f81
--- /dev/null
+++ b/examples/point_120cells_nestml/create_inputs.py
@@ -0,0 +1,25 @@
+import numpy as np
+
+from bmtk.utils.reports.spike_trains import SpikeTrains, PoissonSpikeGenerator
+
+
+def create_inputs_const(pop, firing_rate_hz=10.0, times=(0.0, 3.0)):
+    psg = PoissonSpikeGenerator()  # Uses 'seed' to ensure same results every time
+    psg.add(node_ids='network/{}_nodes.h5'.format(pop), firing_rate=firing_rate_hz, times=times, population=pop)
+    psg.to_sonata('inputs/{}_spikes.h5'.format(pop))
+    # psg.to_csv('inputs/{}_spikes.csv'.format(pop))
+
+
+def create_inputs_sigmoid(pop, min_rate=1.0, max_rate=10.0, onset=0.5, ts_begin=0.0, ts_end=3.0):
+    times = np.linspace(ts_begin, ts_end, 1000)
+    rates = (max_rate-min_rate)/(1.0 + np.exp(-(times-onset)*20)) + min_rate
+
+    psg = PoissonSpikeGenerator()
+    psg.add(node_ids='network/{}_nodes.h5'.format(pop), firing_rate=rates, times=times, population=pop)
+    psg.to_sonata('inputs/{}_spikes.h5'.format(pop))
+    # psg.to_csv('inputs/{}_spikes.csv'.format(pop)
+
+
+if __name__ == '__main__':
+    create_inputs_const(pop='thalamus', firing_rate_hz=10.0)
+    # create_inputs_sigmoid()
diff --git a/examples/point_120cells_nestml/inputs/thalamus_spikes.h5 b/examples/point_120cells_nestml/inputs/thalamus_spikes.h5
new file mode 100644
index 000000000..48cc4d498
Binary files /dev/null and b/examples/point_120cells_nestml/inputs/thalamus_spikes.h5 differ
diff --git a/examples/point_120cells_nestml/network/cortex_node_types.csv b/examples/point_120cells_nestml/network/cortex_node_types.csv
new file mode 100644
index 000000000..1604abf3d
--- /dev/null
+++ b/examples/point_120cells_nestml/network/cortex_node_types.csv
@@ -0,0 +1,3 @@
+node_type_id ei dynamics_params model_type model_template pop_name
+100 e iaf_psc_delta_exc.json point_neuron nestml:iaf_psc_delta LIF_exc
+101 i iaf_psc_delta_inh.json point_neuron nestml:iaf_psc_delta LIF_inh
diff --git a/examples/point_120cells_nestml/plot_output.py b/examples/point_120cells_nestml/plot_output.py
new file mode 100644
index 000000000..bcac1e0d7
--- /dev/null
+++ b/examples/point_120cells_nestml/plot_output.py
@@ -0,0 +1,15 @@
+import matplotlib.pyplot as plt
+
+from bmtk.analyzer.compartment import plot_traces
+from bmtk.analyzer.spike_trains import plot_raster, plot_rates_boxplot
+
+
+# Setting show to False so we can display all the plots at the same time
+plot_raster(config_file='config.simulation.json', group_by='pop_name', show=False)
+plot_rates_boxplot(config_file='config.simulation.json', group_by='pop_name', show=False)
+plot_traces(
+    config_file='config.simulation.json', report_name='membrane_potential', group_by='pop_name',
+    times=(0.0, 200.0), show=False
+)
+
+plt.show()
diff --git a/examples/point_120cells_nestml/run_pointnet.py b/examples/point_120cells_nestml/run_pointnet.py
new file mode 100644
index 000000000..16eaa7ae0
--- /dev/null
+++ b/examples/point_120cells_nestml/run_pointnet.py
@@ -0,0 +1,18 @@
+import os, sys
+
+from bmtk.simulator import pointnet
+
+
+def run(config_file):
+    configure = pointnet.Config.from_json(config_file)
+    configure.build_env()
+
+    network = pointnet.PointNetwork.from_config(configure)
+    sim = pointnet.PointSimulator.from_config(configure, network)
+    sim.run()
+
+
+if __name__ == '__main__':
+    # Find the appropriate config.json file
+    run('config.simulation.json')
+    # run('config.simulation_perturbations.json')
diff --git a/examples/point_components/cell_models/iaf_psc_alpha.nestml b/examples/point_components/cell_models/iaf_psc_alpha.nestml
new file mode 100644
index 000000000..af90787dd
--- /dev/null
+++ b/examples/point_components/cell_models/iaf_psc_alpha.nestml
@@ -0,0 +1,119 @@
+"""
+iaf_psc_alpha - Leaky integrate-and-fire neuron model
+#####################################################
+
+Description
++++++++++++
+
+iaf_psc_alpha is an implementation of a leaky integrate-and-fire model
+with alpha-function kernel synaptic currents. Thus, synaptic currents
+and the resulting post-synaptic potentials have a finite rise time.
+
+The threshold crossing is followed by an absolute refractory period
+during which the membrane potential is clamped to the resting potential.
+
+The general framework for the consistent formulation of systems with
+neuron like dynamics interacting by point events is described in
+[1]_.  A flow chart can be found in [2]_.
+
+Critical tests for the formulation of the neuron model are the
+comparisons of simulation results for different computation step
+sizes.
+
+The iaf_psc_alpha is the standard model used to check the consistency
+of the nest simulation kernel because it is at the same time complex
+enough to exhibit non-trivial dynamics and simple enough compute
+relevant measures analytically.
+
+.. note::
+   If tau_m is very close to tau_syn_exc or tau_syn_inh, numerical problems
+   may arise due to singularities in the propagator matrics. If this is
+   the case, replace equal-valued parameters by a single parameter.
+
+   For details, please see ``IAF_neurons_singularity.ipynb`` in
+   the NEST source code (``docs/model_details``).
+
+
+References
+++++++++++
+
+.. [1] Rotter S,  Diesmann M (1999). Exact simulation of
+       time-invariant linear systems with applications to neuronal
+       modeling. Biologial Cybernetics 81:381-402.
+       DOI: https://doi.org/10.1007/s004220050570
+.. [2] Diesmann M, Gewaltig M-O, Rotter S, & Aertsen A (2001). State
+       space analysis of synchronous spiking in cortical neural
+       networks. Neurocomputing 38-40:565-571.
+       DOI: https://doi.org/10.1016/S0925-2312(01)00409-X
+.. [3] Morrison A, Straube S, Plesser H E, Diesmann M (2006). Exact
+       subthreshold integration with continuous spike times in discrete time
+       neural network simulations. Neural Computation, in press
+       DOI: https://doi.org/10.1162/neco.2007.19.1.47
+
+
+See also
+++++++++
+
+iaf_psc_delta, iaf_psc_exp, iaf_cond_alpha
+"""
+neuron iaf_psc_alpha:
+
+  state:
+    r integer = 0                 # counts number of tick during the refractory period
+    V_abs mV = 0 mV
+  end
+
+  equations:
+    kernel I_kernel_inh = (e / tau_syn_inh) * t * exp(-t / tau_syn_inh)
+    kernel I_kernel_exc = (e / tau_syn_exc) * t * exp(-t / tau_syn_exc)
+    recordable inline V_m mV = V_abs + E_L    # Membrane potential
+    inline I pA = convolve(I_kernel_exc, exc_spikes) - convolve(I_kernel_inh, inh_spikes) + I_e + I_stim
+    V_abs' = -V_abs / tau_m + I / C_m
+  end
+
+  parameters:
+    C_m     pF = 250 pF       # Capacitance of the membrane
+    tau_m   ms = 10 ms        # Membrane time constant
+    tau_syn_inh ms = 2 ms     # Time constant of synaptic current
+    tau_syn_exc ms = 2 ms     # Time constant of synaptic current
+    t_ref   ms = 2 ms         # Duration of refractory period
+    E_L     mV = -70 mV       # Resting potential
+    V_reset mV = -70 mV - E_L # Reset potential of the membrane
+    V_th    mV = -55 mV - E_L  # Spike threshold
+
+    # constant external input current
+    I_e pA = 0 pA
+  end
+
+  internals:
+    RefractoryCounts integer = steps(t_ref)    # refractory time in steps
+  end
+
+  input:
+    exc_spikes pA <- excitatory spike
+    inh_spikes pA <- inhibitory spike
+    I_stim pA <- continuous
+  end
+
+  output: spike
+
+  update:
+    if r == 0: # neuron not refractory
+      integrate_odes()
+    else: # neuron is absolute refractory
+      r = r - 1
+    end
+
+    if V_abs >= V_th: # threshold crossing
+      # A supra-threshold membrane potential should never be observable.
+      # The reset at the time of threshold crossing enables accurate
+      # integration independent of the computation step size, see [2,3] for
+      # details.
+      r = RefractoryCounts
+      V_abs = V_reset
+      emit_spike()
+    end
+
+  end
+
+end
diff --git a/examples/point_components/cell_models/iaf_psc_delta.nestml b/examples/point_components/cell_models/iaf_psc_delta.nestml
new file mode 100644
index 000000000..58c60ac79
--- /dev/null
+++ b/examples/point_components/cell_models/iaf_psc_delta.nestml
@@ -0,0 +1,125 @@
+"""
+iaf_psc_delta - Current-based leaky integrate-and-fire neuron model with delta-kernel post-synaptic currents
+############################################################################################################
+
+Description
++++++++++++
+
+iaf_psc_delta is an implementation of a leaky integrate-and-fire model
+where the potential jumps on each spike arrival.
+
+The threshold crossing is followed by an absolute refractory period
+during which the membrane potential is clamped to the resting potential.
+
+Spikes arriving while the neuron is refractory, are discarded by
+default. If the property ``with_refr_input`` is set to true, such
+spikes are added to the membrane potential at the end of the
+refractory period, dampened according to the interval between
+arrival and end of refractoriness.
+
+The general framework for the consistent formulation of systems with
+neuron like dynamics interacting by point events is described in
+[1]_.  A flow chart can be found in [2]_.
+
+Critical tests for the formulation of the neuron model are the
+comparisons of simulation results for different computation step
+sizes. sli/testsuite/nest contains a number of such tests.
+
+The iaf_psc_delta is the standard model used to check the consistency
+of the nest simulation kernel because it is at the same time complex
+enough to exhibit non-trivial dynamics and simple enough compute
+relevant measures analytically.
+
+
+References
+++++++++++
+
+.. [1] Rotter S,  Diesmann M (1999). Exact simulation of
+       time-invariant linear systems with applications to neuronal
+       modeling. Biologial Cybernetics 81:381-402.
+       DOI: https://doi.org/10.1007/s004220050570
+.. [2] Diesmann M, Gewaltig M-O, Rotter S, & Aertsen A (2001). State
+       space analysis of synchronous spiking in cortical neural
+       networks. Neurocomputing 38-40:565-571.
+       DOI: https://doi.org/10.1016/S0925-2312(01)00409-X
+
+
+See also
+++++++++
+
+iaf_psc_alpha, iaf_psc_exp
+"""
+neuron iaf_psc_delta:
+
+  state:
+    refr_spikes_buffer mV = 0 mV
+    r integer = 0 # counts number of tick during the refractory period
+    V_abs mV = 0 mV
+  end
+
+  equations:
+    kernel G = delta(t)
+    recordable inline V_m mV = V_abs + E_L # Membrane potential.
+    V_abs' = -V_abs / tau_m + (mV / pA / ms) * convolve(G, spikes) + (I_e + I_stim) / C_m
+  end
+
+  parameters:
+    tau_m   ms = 10 ms      # Membrane time constant.
+    C_m     pF = 250 pF     # Capacity of the membrane
+    t_ref   ms = 2 ms       # Duration of refractory period.
+    tau_syn ms = 2 ms       # Time constant of synaptic current.
+    E_L     mV = -70 mV     # Resting membrane potential.
+    V_reset mV = -70 mV - E_L # Reset potential of the membrane.
+    Theta   mV = -55 mV - E_L # Spike threshold.
+    V_min mV = -inf * 1 mV           # Absolute lower value for the membrane potential
+    with_refr_input boolean = false # If true, do not discard input during  refractory period. Default: false.
+
+    # constant external input current
+    I_e pA = 0 pA
+  end
+
+  internals:
+    h ms = resolution()
+    RefractoryCounts integer = steps(t_ref) # refractory time in steps
+  end
+
+  input:
+    spikes pA <- spike
+    I_stim pA <- continuous
+  end
+
+  output: spike
+
+  update:
+    if r == 0: # neuron not refractory
+      integrate_odes()
+
+      # if we have accumulated spikes from refractory period,
+      # add and reset accumulator
+      if with_refr_input and refr_spikes_buffer != 0.0 mV:
+        V_abs += refr_spikes_buffer
+        refr_spikes_buffer = 0.0 mV
+      end
+
+      # lower bound of membrane potential
+      V_abs = V_abs < V_min?V_min:V_abs
+
+    else: # neuron is absolute refractory
+      # read spikes from buffer and accumulate them, discounting
+      # for decay until end of refractory period
+      # the buffer is clear automatically
+      if with_refr_input:
+        refr_spikes_buffer += spikes * exp(-r * h / tau_m) * mV/pA
+      end
+      r -= 1
+    end
+
+    if V_abs >= Theta: # threshold crossing
+        r = RefractoryCounts
+        V_abs = V_reset
+        emit_spike()
+    end
+
+  end
+
+end
diff --git a/examples/point_components/synaptic_models/static_synapse.nestml b/examples/point_components/synaptic_models/static_synapse.nestml
new file mode 100644
index 000000000..8a30b14aa
--- /dev/null
+++ b/examples/point_components/synaptic_models/static_synapse.nestml
@@ -0,0 +1,26 @@
+"""
+Static synapse
+##############
+
+Description
++++++++++++
+A synapse where the synaptic strength (weight) does not evolve with simulated time, but is defined as a (constant) parameter.
+"""
+synapse static:
+
+  parameters:
+    w real = 900  @nest::weight @homogeneous
+    d ms = .9 ms  @nest::delay @heterogeneous
+    a real = 3.141592653589793  @nest::a @homogeneous
+    b real = 100.  @nest::b @heterogeneous
+  end
+
+  input:
+    pre_spikes mV <- spike
+  end
+
+  onReceive(pre_spikes):
+    deliver_spike(3.18E-3 * a * b * w, d)
+  end
+
+end