diff --git a/pyvcell/data_model/simulation.py b/pyvcell/data_model/simulation.py
index 5a15830..a54f27e 100644
--- a/pyvcell/data_model/simulation.py
+++ b/pyvcell/data_model/simulation.py
@@ -7,6 +7,7 @@
 from pyvcell.api.vcell_client import ApiClient, ApiResponse, Configuration, SolverResourceApi
 from pyvcell.data_model.result import Result
 from pyvcell.data_model.spatial_model import SpatialModel
+from pyvcell.data_model.vcml_spatial_model import VcmlSpatialModel
 from pyvcell.solvers.fvsolver import solve as fvsolve
 
 
@@ -40,7 +41,9 @@ def run(self, duration: float | None = None, output_time_step: float | None = No
         # create temp file to write sbml document to
         sbml_path = self.out_dir / "model.xml"
         self.model.export(sbml_path)
-        response: ApiResponse[bytearray] = solver_api.get_fv_solver_input_from_sbml_with_http_info(str(sbml_path))
+        response: ApiResponse[bytearray] = solver_api.get_fv_solver_input_from_sbml_with_http_info(
+            str(sbml_path), duration=duration, output_time_step=output_time_step
+        )
         sbml_path.unlink()
         if response.status_code != 200:
             raise ValueError(f"Failed to get solver input files: {response}")
@@ -69,3 +72,60 @@ def run(self, duration: float | None = None, output_time_step: float | None = No
 
     def cleanup(self) -> None:
         shutil.rmtree(self.out_dir)
+
+
+class VcmlSpatialSimulation:
+    model: VcmlSpatialModel
+    out_dir: Path
+
+    def __init__(self, model: VcmlSpatialModel, out_dir: Path | None = None):
+        self.model = model
+        if out_dir is None:
+            self.out_dir = Path(tempfile.mkdtemp(prefix="out_dir_"))
+        else:
+            self.out_dir = out_dir
+
+    def run(self, simulation_name: str) -> Result:
+        # create an unauthenticated API client
+        api_url: str = "https://vcell-dev.cam.uchc.edu"  # vcell base url
+        api_client = ApiClient(Configuration(host=api_url))
+        solver_api = SolverResourceApi(api_client)
+
+        # prepare solver input files
+        # 1. upload the SBML model and retrieve generated solver inputs as a zip file
+        # 2. extract the zip archive into the output directory
+        # 3. remove the zip archive
+        # create temp file to write sbml document to
+        vcml_path = self.out_dir / "model.xml"
+        self.model.export(vcml_path)
+        response: ApiResponse[bytearray] = solver_api.get_fv_solver_input_from_vcml_with_http_info(
+            vcml_file=str(vcml_path), simulation_name=simulation_name
+        )
+        vcml_path.unlink()
+        if response.status_code != 200:
+            raise ValueError(f"Failed to get solver input files: {response}")
+        zip_archive = self.out_dir / "solver_input_files.zip"
+        with open(zip_archive, "wb") as f:
+            f.write(response.data)
+        shutil.unpack_archive(zip_archive, self.out_dir)
+        zip_archive.unlink()
+
+        # identify sim_id and job_id from the solver input files
+        files: list[str] = os.listdir(self.out_dir)
+        fv_input_file: Path | None = next((self.out_dir / file for file in files if file.endswith(".fvinput")), None)
+        vcg_input_file: Path | None = next((self.out_dir / file for file in files if file.endswith(".vcg")), None)
+        if fv_input_file is None or vcg_input_file is None:
+            raise ValueError(".fvinput file or .vcg file not found")
+        sim_id = int(fv_input_file.name.split("_")[1])
+        job_id = int(fv_input_file.name.split("_")[2])
+
+        # run the simulation
+        ret_code = fvsolve(input_file=fv_input_file, vcg_file=vcg_input_file, output_dir=self.out_dir)
+        if ret_code != 0:
+            raise ValueError(f"Error in solve: {ret_code}")
+
+        # return the result
+        return Result(solver_output_dir=self.out_dir, sim_id=sim_id, job_id=job_id)
+
+    def cleanup(self) -> None:
+        shutil.rmtree(self.out_dir)
diff --git a/pyvcell/data_model/vcml_spatial_model.py b/pyvcell/data_model/vcml_spatial_model.py
new file mode 100644
index 0000000..9b442f1
--- /dev/null
+++ b/pyvcell/data_model/vcml_spatial_model.py
@@ -0,0 +1,64 @@
+import os
+from pathlib import Path
+from typing import Any, Union
+
+import pyvcell.vcml as vc
+
+
+class VcmlSpatialModel:
+    """
+    Spatial extension of `libsbml.Model`. All class methods are inherited from `libsbml.Model`: see libsbml documentation for more details.
+    This class is constructed with one of 3 entrypoints: either the filepath to a valid SBMLSpatial model, OR level, version, model_id, OR model_id
+    """
+
+    _bio_model: vc.Biomodel
+
+    def __init__(self, filepath: Path) -> None:
+        reader: vc.VcmlReader = vc.VcmlReader()
+        # read filepath as string
+        with open(filepath, "r") as file:
+            self._bio_model = reader.parse_biomodel(file.read())
+
+    @property
+    def bio_model(self) -> vc.Biomodel:
+        return self._bio_model
+
+    @property
+    def model(self) -> vc.Model:
+        if self._bio_model.model is None:
+            raise ValueError("Model is not defined in the VCML document.")
+        return self._bio_model.model
+
+    def export(self, filename: Union[os.PathLike[str], str]) -> None:
+        document = vc.VCMLDocument(biomodel=self._bio_model)
+        vcml_str: str = vc.VcmlWriter().write_vcml(document=document)
+        with open(filename, "w") as file:
+            file.write(vcml_str)
+
+    @property
+    def model_parameters(self) -> list[vc.ModelParameter]:
+        if self._bio_model.model is None:
+            raise ValueError("Model is not defined in the VCML document.")
+        return self._bio_model.model.model_parameters
+
+    @property
+    def kinetics_parameters(self) -> list[vc.KineticsParameter]:
+        kinetics_parameters: list[vc.KineticsParameter] = []
+        for reaction in self.model.reactions:
+            if reaction.kinetics is None:
+                continue
+            for param in reaction.kinetics.kinetics_parameters:
+                param.reaction_name = reaction.name
+                kinetics_parameters.append(param)
+        return kinetics_parameters
+
+    @property
+    def applications(self) -> list[vc.Application]:
+        return self.bio_model.applications
+
+    @property
+    def simulations(self) -> list[vc.Simulation]:
+        sims: list[vc.Simulation] = []
+        for app in self.applications:
+            sims.extend(app.simulations)
+        return sims
diff --git a/pyvcell/vcml/__init__.py b/pyvcell/vcml/__init__.py
index 9fdbbb9..1d3cf82 100644
--- a/pyvcell/vcml/__init__.py
+++ b/pyvcell/vcml/__init__.py
@@ -9,6 +9,7 @@
     Model,
     ModelParameter,
     Reaction,
+    Simulation,
     Species,
     SpeciesMapping,
     SpeciesReference,
@@ -40,4 +41,6 @@
     "SurfaceClass",
     "SpeciesMapping",
     "BoundaryType",
+    "Application",
+    "Simulation",
 ]
diff --git a/pyvcell/vcml/models.py b/pyvcell/vcml/models.py
index f348109..19025c9 100644
--- a/pyvcell/vcml/models.py
+++ b/pyvcell/vcml/models.py
@@ -33,7 +33,7 @@ class ModelParameter(Parameter):
 
 
 class KineticsParameter(Parameter):
-    pass
+    reaction_name: str
 
 
 class Kinetics(VcmlNode):
@@ -145,6 +145,13 @@ class ReactionMapping(VcmlNode):
     included: bool = True
 
 
+class Simulation(VcmlNode):
+    name: str
+    duration: float
+    output_time_step: float
+    mesh_size: tuple[int, int, int]
+
+
 class Application(VcmlNode):
     name: str
     stochastic: bool
@@ -152,6 +159,7 @@ class Application(VcmlNode):
     compartment_mappings: list[CompartmentMapping] = Field(default_factory=list)
     species_mappings: list[SpeciesMapping] = Field(default_factory=list)
     reaction_mappings: list[ReactionMapping] = Field(default_factory=list)
+    simulations: list[Simulation] = Field(default_factory=list)
 
 
 class Biomodel(VcmlNode):
diff --git a/pyvcell/vcml/vcml_reader.py b/pyvcell/vcml/vcml_reader.py
index b61096e..07c83a7 100644
--- a/pyvcell/vcml/vcml_reader.py
+++ b/pyvcell/vcml/vcml_reader.py
@@ -26,13 +26,15 @@ def strip_namespace(tag: str) -> str:
 
 class VcmlReader:
     @staticmethod
-    def parse_biomodel(xml_string: str) -> vc.Biomodel | None:
+    def parse_biomodel(xml_string: str) -> vc.Biomodel:
         xml_string = xml_string.replace('<?xml version="1.0" encoding="UTF-8"?>', "")
         xml_string = xml_string.replace("<?xml version='1.0' encoding='UTF-8'?>", "")
         root = etree.fromstring(xml_string)
         document = vc.VCMLDocument()
         visitor = BiomodelVisitor(document)
         visitor.visit(root, document)
+        if visitor.document.biomodel is None:
+            raise ValueError("No biomodel found")
         return visitor.document.biomodel
 
     @classmethod
@@ -135,7 +137,13 @@ def visit_Parameter(self, element: _Element, node: vc.Model | vc.Kinetics) -> No
             parameter = model_parameter
         elif strip_namespace(parent.tag) == "Kinetics":
             kinetics: vc.Kinetics = node  # type: ignore[assignment]
-            kinetics_parameter = vc.KineticsParameter(name=name, value=value, role=role, unit=unit)
+            reaction_node = parent.getparent()
+            if reaction_node is None:
+                raise ValueError("Kinetics element has no parent")
+            reaction_name = reaction_node.get("Name", default="unknown")
+            kinetics_parameter = vc.KineticsParameter(
+                name=name, value=value, role=role, unit=unit, reaction_name=reaction_name
+            )
             kinetics.kinetics_parameters.append(kinetics_parameter)
             parameter = kinetics_parameter
         else:
@@ -150,6 +158,32 @@ def visit_SimulationSpec(self, element: _Element, node: vc.Biomodel) -> None:
         node.applications.append(application)
         self.generic_visit(element, application)
 
+    def visit_Simulation(self, element: _Element, node: vc.Application) -> None:
+        name: str = element.get("Name", default="unnamed")
+        duration: float | None = None
+        output_time_step: float | None = None
+        mesh_size: tuple[int, int, int] | None = None
+        for sim_child in element:
+            if strip_namespace(sim_child.tag) == "SolverTaskDescription":
+                solver_task_description_element = sim_child
+                for child in solver_task_description_element:
+                    if strip_namespace(child.tag) == "TimeBound":
+                        duration = float(child.get("EndTime", default="5.0"))
+                    elif strip_namespace(child.tag) == "OutputOptions":
+                        output_time_step = float(child.get("OutputTimeStep", default="0.1"))
+            elif strip_namespace(sim_child.tag) == "MeshSpecification":
+                mesh_specification_element = sim_child
+                for mesh_child in mesh_specification_element:
+                    if strip_namespace(mesh_child.tag) == "Size":
+                        mesh_x = int(mesh_child.get("X", default="1"))
+                        mesh_y = int(mesh_child.get("Y", default="1"))
+                        mesh_z = int(mesh_child.get("Z", default="1"))
+                        mesh_size = (mesh_x, mesh_y, mesh_z)
+        if duration is None or output_time_step is None or mesh_size is None:
+            raise ValueError("Simulation element is missing required child elements")
+        simulation = vc.Simulation(name=name, duration=duration, output_time_step=output_time_step, mesh_size=mesh_size)
+        node.simulations.append(simulation)
+
     def visit_Geometry(self, element: _Element, node: vc.Application) -> None:
         name: str = element.get("Name", default="unnamed")
         dim = int(element.get("Dimension", default="0"))
diff --git a/pyvcell/vcml/vcml_writer.py b/pyvcell/vcml/vcml_writer.py
index 53ad3ce..fbb7d52 100644
--- a/pyvcell/vcml/vcml_writer.py
+++ b/pyvcell/vcml/vcml_writer.py
@@ -36,11 +36,19 @@ def write_model(self, model: vc.Model, parent: _Element) -> None:
             parameter_element = Element("Parameter", Name=parameter.name, Role=parameter.role, Unit=parameter.unit)
             parameter_element.text = str(parameter.value)
             model_parameters_element.append(parameter_element)
+        for species in model.species:
+            species_type_element = Element("Compound", Name=species.name)
+            annotation_element = Element("Annotation")
+            annotation_element.text = species.name
+            species_type_element.append(annotation_element)
+            parent.append(species_type_element)
         for compartment in model.compartments:
             compartment_element = Element("Feature" if compartment.dim == 3 else "Membrane", Name=compartment.name)
             parent.append(compartment_element)
         for species in model.species:
-            species_element = Element("LocalizedCompound", Name=species.name, Structure=species.structure_name)
+            species_element = Element(
+                "LocalizedCompound", Name=species.name, CompoundRef=species.name, Structure=species.structure_name
+            )
             parent.append(species_element)
         for reaction in model.reactions:
             reaction_element = Element("SimpleReaction", Structure=reaction.compartment_name, Name=reaction.name)
@@ -73,6 +81,8 @@ def write_application(self, application: vc.Application, parent: _Element) -> No
         geometry_element = Element("Geometry", Name=application.geometry.name, Dimension=str(application.geometry.dim))
         parent.append(geometry_element)
         self.write_geometry(application.geometry, geometry_element)
+
+        # ---- geometry context -----
         geometry_context_element = Element("GeometryContext")
         parent.append(geometry_context_element)
         for compartment_mapping in application.compartment_mappings:
@@ -94,9 +104,13 @@ def write_application(self, application: vc.Application, parent: _Element) -> No
             )
             mapping_element.append(boundaries_types_element)
             geometry_context_element.append(mapping_element)
+
+        # ---- reaction context -----
+        reaction_context_element = Element("ReactionContext")
+        parent.append(reaction_context_element)
         for species_mapping in application.species_mappings:
             mapping_element = Element("LocalizedCompoundSpec", LocalizedCompoundRef=species_mapping.species_name)
-            parent.append(mapping_element)
+            reaction_context_element.append(mapping_element)
             self.write_species_mapping(species_mapping, mapping_element)
         for reaction_mapping in application.reaction_mappings:
             mapping_element = Element(
@@ -104,7 +118,50 @@ def write_application(self, application: vc.Application, parent: _Element) -> No
                 ReactionStepRef=reaction_mapping.reaction_name,
                 ReactionMapping="included" if reaction_mapping.included else "excluded",
             )
-            parent.append(mapping_element)
+            reaction_context_element.append(mapping_element)
+
+        # ---- mathDescription ---- (skip this for now)
+        math_description_element = Element("MathDescription", Name="dummy_math_description")
+        parent.append(math_description_element)
+
+        # ---- simulations -----
+        for simulation in application.simulations:
+            simulation_element = Element("Simulation", Name=simulation.name)
+            parent.append(simulation_element)
+            solver_task_description_element = Element(
+                "SolverTaskDescription",
+                TaskType="Unsteady",
+                UseSymbolicJacobian="false",
+                Solver="Sundials Stiff PDE Solver (Variable Time Step)",
+            )
+            simulation_element.append(solver_task_description_element)
+            solver_task_description_element.append(
+                Element("TimeBound", StartTime="0.0", EndTime=str(simulation.duration))
+            )
+            solver_task_description_element.append(
+                Element("TimeStep", DefaultTime="0.05", MinTime="0.0", MaxTime="0.1")
+            )
+            solver_task_description_element.append(Element("ErrorTolerance", Absolut="1.0E-9", Relative="1.0E-7"))
+            solver_task_description_element.append(
+                Element("OutputOptions", OutputTimeStep=str(simulation.output_time_step))
+            )
+
+            sundials_solver_options_element = Element("SundialsSolverOptions")
+            max_order_advection_element = Element("maxOrderAdvection")
+            max_order_advection_element.text = "2"
+            sundials_solver_options_element.append(max_order_advection_element)
+            solver_task_description_element.append(sundials_solver_options_element)
+            number_processors_element = Element("NumberProcessors")
+            number_processors_element.text = "1"
+            solver_task_description_element.append(number_processors_element)
+            simulation_element.append(Element("MathOverrides"))
+
+            mesh_specification_element = Element("MeshSpecification")
+            size_element = Element(
+                "Size", X=str(simulation.mesh_size[0]), Y=str(simulation.mesh_size[1]), Z=str(simulation.mesh_size[2])
+            )
+            mesh_specification_element.append(size_element)
+            simulation_element.append(mesh_specification_element)
 
     def write_geometry(self, geometry: vc.Geometry, parent: _Element) -> None:
         extent_element = Element(
diff --git a/tests/conftest.py b/tests/conftest.py
index 19bc9da..49efcc5 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,3 +1,3 @@
 from tests.fixtures.data_fixtures import solver_output_path, solver_output_simid_jobid, zarr_path  # noqa: F401
 from tests.fixtures.model_fixtures import sbml_spatial_model_1d_path, sbml_spatial_model_3d_path  # noqa: F401
-from tests.fixtures.vcell_model_fixtures import vcml_spatial_model_1d_path  # noqa: F401
+from tests.fixtures.vcell_model_fixtures import vcml_spatial_model_1d_path, vcml_spatial_small_3d_path  # noqa: F401
diff --git a/tests/data_model/test_simulations.py b/tests/data_model/test_simulations.py
index 1e33912..c20452c 100644
--- a/tests/data_model/test_simulations.py
+++ b/tests/data_model/test_simulations.py
@@ -2,12 +2,14 @@
 
 import numpy as np
 
+import pyvcell.vcml as vc
 from pyvcell.data_model.result import Result
-from pyvcell.data_model.simulation import SpatialSimulation
+from pyvcell.data_model.simulation import SpatialSimulation, VcmlSpatialSimulation
 from pyvcell.data_model.spatial_model import SpatialModel
+from pyvcell.data_model.vcml_spatial_model import VcmlSpatialModel
 
 
-def test_model_parse_1d(sbml_spatial_model_1d_path: Path) -> None:
+def test_sbml_model_parse_1d(sbml_spatial_model_1d_path: Path) -> None:
     assert sbml_spatial_model_1d_path.is_file()
 
     spatial_model = SpatialModel(filepath=sbml_spatial_model_1d_path)
@@ -29,14 +31,16 @@ def test_model_parse_1d(sbml_spatial_model_1d_path: Path) -> None:
     assert spatial_model.get_coordinate_symbols() == ["x"]
 
     simulation_orig = SpatialSimulation(model=spatial_model)
-    results_orig: Result = simulation_orig.run(duration=5.0, output_time_step=0.1)
+    # results_orig: Result = simulation_orig.run(duration=5.0, output_time_step=0.1)
+    results_orig: Result = simulation_orig.run()
 
     spatial_model.set_parameter_value("Kr_r0", 100.0)
     assert spatial_model.model.getParameter("Kr_r0").getValue() == 100.0
     spatial_model.set_parameter_value("s0_BC_Xm", 1.0)
 
     simulation_changed = SpatialSimulation(model=spatial_model)
-    results_changed: Result = simulation_changed.run(duration=5.0, output_time_step=0.1)
+    #    results_changed: Result = simulation_changed.run(duration=5.0, output_time_step=0.1)
+    results_changed: Result = simulation_changed.run()
 
     channels_orig = results_orig.channels
     channels_changed = results_changed.channels
@@ -65,7 +69,7 @@ def test_model_parse_1d(sbml_spatial_model_1d_path: Path) -> None:
     )
 
 
-def test_model_parse_3d(sbml_spatial_model_3d_path: Path) -> None:
+def test_sbml_model_parse_3d(sbml_spatial_model_3d_path: Path) -> None:
     assert sbml_spatial_model_3d_path.is_file()
 
     spatial_model = SpatialModel(filepath=sbml_spatial_model_3d_path)
@@ -114,13 +118,15 @@ def test_model_parse_3d(sbml_spatial_model_3d_path: Path) -> None:
     }
 
     simulation_orig = SpatialSimulation(model=spatial_model)
-    results_orig: Result = simulation_orig.run(duration=5.0, output_time_step=0.1)
+    # results_orig: Result = simulation_orig.run(duration=5.0, output_time_step=0.1)
+    results_orig: Result = simulation_orig.run()
 
     spatial_model.set_parameter_value("Kr_r0", 100.0)
     assert spatial_model.model.getParameter("Kr_r0").getValue() == 100.0
 
     simulation_changed = SpatialSimulation(model=spatial_model)
-    results_changed: Result = simulation_changed.run(duration=5.0, output_time_step=0.1)
+    # results_changed: Result = simulation_changed.run(duration=5.0, output_time_step=0.1)
+    results_changed: Result = simulation_changed.run()
 
     channels_orig = results_orig.channels
     channels_changed = results_changed.channels
@@ -170,3 +176,44 @@ def test_model_parse_3d(sbml_spatial_model_3d_path: Path) -> None:
             dtype=np.float64,
         ),
     )
+
+
+def test_vcml_model_parse_3d(vcml_spatial_model_1d_path: Path) -> None:
+    assert vcml_spatial_model_1d_path.is_file()
+
+    vcml_spatial_model = VcmlSpatialModel(filepath=vcml_spatial_model_1d_path)
+    assert vcml_spatial_model is not None
+    parameters: list[vc.ModelParameter] = vcml_spatial_model.model_parameters
+    assert {p.name: p.value for p in parameters} == {"Kf_r0": 1.0, "Kr_r0": 0.5}
+
+    sim_name = vcml_spatial_model.applications[0].simulations[0].name
+    simulation_orig = VcmlSpatialSimulation(model=vcml_spatial_model)
+    # # create a temporary file to write the VCML content to
+    # vcml_path = Path(tempfile.mktemp(prefix="model_", suffix=".xml"))
+    # vcml_spatial_model.export(vcml_path)
+    results_orig: Result = simulation_orig.run(simulation_name=sim_name)
+
+    Kr_r0: vc.ModelParameter = next(p for p in vcml_spatial_model.model_parameters if p.name == "Kr_r0")
+    Kr_r0.value = 100.0
+
+    simulation_changed = VcmlSpatialSimulation(model=vcml_spatial_model)
+    results_changed: Result = simulation_changed.run(simulation_name=sim_name)
+
+    channels_orig = results_orig.channels
+    channels_changed = results_changed.channels
+    assert [channel.label for channel in channels_orig] == ["s0", "s1", "J_r0", "s0_init_uM", "s1_init_uM"]
+    assert [channel.label for channel in channels_changed] == ["s0", "s1", "J_r0", "s0_init_uM", "s1_init_uM"]
+    assert np.allclose(
+        results_orig.concentrations[0, 0::10],
+        np.array(
+            [500000.00000000006, 236185.00811512678, 111567.6111715472],
+            dtype=np.float64,
+        ),
+    )
+    assert np.allclose(
+        results_changed.concentrations[0, 0::10],
+        np.array(
+            [500000.00000000006, 9.900990099017005, 9.900990099052574],
+            dtype=np.float64,
+        ),
+    )
diff --git a/tests/fixtures/data/SmallSpacialProject_3D.vcml b/tests/fixtures/data/SmallSpacialProject_3D.vcml
new file mode 100644
index 0000000..0048bcc
--- /dev/null
+++ b/tests/fixtures/data/SmallSpacialProject_3D.vcml
@@ -0,0 +1,184 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--This biomodel was generated in VCML Version Alpha_Version_7.7.0_build_15-->
+<vcml xmlns="http://sourceforge.net/projects/vcell/vcml" Version="Alpha_Version_7.7.0_build_15">
+  <BioModel Name="TinySpacialProject_Application0_unnamed_spatialGeom">
+    <Model Name="unnamed">
+      <ModelParameters>
+        <Parameter Name="Kf_r0" Role="user defined" Unit="s-1" SbmlName="Kf_r0">1.0</Parameter>
+        <Parameter Name="Kr_r0" Role="user defined" Unit="s-1" SbmlName="Kr_r0">0.5</Parameter>
+        <Parameter Name="Kf_r1" Role="user defined" Unit="l.dm-2.s-1">1.0</Parameter>
+        <Parameter Name="Kr_r1" Role="user defined" Unit="s-1">1.0</Parameter>
+        <Parameter Name="Kf_r2" Role="user defined" Unit="s-1">1.0</Parameter>
+        <Parameter Name="Kr_r2" Role="user defined" Unit="l.dm-2.s-1">1.0</Parameter>
+      </ModelParameters>
+      <Compound Name="s0">
+        <Annotation>s0</Annotation>
+      </Compound>
+      <Compound Name="s1">
+        <Annotation>s1</Annotation>
+      </Compound>
+      <Compound Name="s2">
+        <Annotation>s2</Annotation>
+      </Compound>
+      <Compound Name="s3">
+        <Annotation>s3</Annotation>
+      </Compound>
+      <Feature Name="c0" />
+      <Feature Name="c1" />
+      <Membrane MembraneVoltage="Voltage_m0" Name="m0" />
+      <LocalizedCompound Name="s0" SbmlName="s0" CompoundRef="s0" Structure="c0" OverrideName="true" />
+      <LocalizedCompound Name="s1" SbmlName="s1" CompoundRef="s1" Structure="c0" OverrideName="true" />
+      <LocalizedCompound Name="s2" SbmlName="s2" CompoundRef="s2" Structure="m0" OverrideName="true" />
+      <LocalizedCompound Name="s3" SbmlName="s3" CompoundRef="s3" Structure="c1" OverrideName="true" />
+      <SimpleReaction Structure="c0" Name="r0" Reversible="true" FluxOption="MolecularOnly" SbmlName="r0">
+        <Reactant LocalizedCompoundRef="s0" Stoichiometry="1" />
+        <Product LocalizedCompoundRef="s1" Stoichiometry="1" />
+        <Kinetics KineticsType="GeneralKinetics">
+          <Parameter Name="J" Role="reaction rate" Unit="umol.l-1.s-1">((Kf_r0 * s0) - (1.0 * Kr_r0 * s1))</Parameter>
+        </Kinetics>
+      </SimpleReaction>
+      <SimpleReaction Structure="m0" Name="r1" Reversible="true" FluxOption="MolecularOnly" SbmlName="r1">
+        <Reactant LocalizedCompoundRef="s0" Stoichiometry="1" />
+        <Product LocalizedCompoundRef="s2" Stoichiometry="1" />
+        <Kinetics KineticsType="GeneralKinetics">
+          <Parameter Name="J" Role="reaction rate" Unit="umol.dm-2.s-1">((Kf_r1 * s0) - (Kr_r1 * s2))</Parameter>
+          <Parameter Name="I" Role="inward current density" Unit="pA.dm-2">0.0</Parameter>
+          <Parameter Name="netValence" Role="net charge valence" Unit="1">1.0</Parameter>
+        </Kinetics>
+      </SimpleReaction>
+      <SimpleReaction Structure="m0" Name="r2" Reversible="true" FluxOption="MolecularOnly" SbmlName="r2">
+        <Reactant LocalizedCompoundRef="s2" Stoichiometry="1" />
+        <Product LocalizedCompoundRef="s3" Stoichiometry="1" />
+        <Kinetics KineticsType="GeneralKinetics">
+          <Parameter Name="J" Role="reaction rate" Unit="umol.dm-2.s-1">((Kf_r2 * s2) - (Kr_r2 * s3))</Parameter>
+          <Parameter Name="I" Role="inward current density" Unit="pA.dm-2">0.0</Parameter>
+          <Parameter Name="netValence" Role="net charge valence" Unit="1">1.0</Parameter>
+        </Kinetics>
+      </SimpleReaction>
+      <Diagram Name="c0" Structure="c0">
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s2" LocationX="47" LocationY="220" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s1" LocationX="14" LocationY="220" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s3" LocationX="14" LocationY="34" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r0" LocationX="111" LocationY="147" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r1" LocationX="14" LocationY="39" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r2" LocationX="74" LocationY="34" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s0" LocationX="170" LocationY="34" />
+      </Diagram>
+      <Diagram Name="m0" Structure="m0">
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s2" LocationX="47" LocationY="220" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s1" LocationX="14" LocationY="220" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s3" LocationX="14" LocationY="34" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r0" LocationX="111" LocationY="147" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r1" LocationX="14" LocationY="39" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r2" LocationX="74" LocationY="34" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s0" LocationX="170" LocationY="34" />
+      </Diagram>
+      <Diagram Name="c1" Structure="c1">
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s2" LocationX="47" LocationY="220" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s1" LocationX="14" LocationY="220" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s3" LocationX="14" LocationY="34" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r0" LocationX="111" LocationY="147" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r1" LocationX="14" LocationY="39" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r2" LocationX="74" LocationY="34" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s0" LocationX="170" LocationY="34" />
+      </Diagram>
+      <ModelUnitSystem VolumeSubstanceUnit="umol" MembraneSubstanceUnit="umol" LumpedReactionSubstanceUnit="umol" VolumeUnit="l" AreaUnit="dm2" LengthUnit="dm" TimeUnit="s" />
+    </Model>
+    <SimulationSpec Name="unnamed_spatialGeom" Stochastic="false" UseConcentration="true" SpringSaLaD="false" RuleBased="false" MassConservationModelReduction="false" InsufficientIterations="false" InsufficientMaxMolecules="false" CharacteristicSize="0.10403160765672961">
+      <NetworkConstraints RbmMaxIteration="1" RbmMaxMoleculesPerSpecies="10" RbmSpeciesLimit="800" RbmReactionsLimit="2500" />
+      <Annotation />
+      <Geometry Name="spatialGeom" Dimension="3">
+        <Extent X="10.0" Y="10.0" Z="10.0" />
+        <Origin X="0.0" Y="0.0" Z="0.0" />
+        <SubVolume Name="subdomain1" Handle="1" Type="Analytical">
+          <AnalyticExpression>((pow((-5.0 + x),2.0) + pow((-5.0 + y),2.0) + pow((-5.0 + z),2.0)) &lt; 16.0)</AnalyticExpression>
+        </SubVolume>
+        <SubVolume Name="subdomain0" Handle="0" Type="Analytical">
+          <AnalyticExpression>1.0</AnalyticExpression>
+        </SubVolume>
+        <SurfaceClass Name="subdomain0_subdomain1_membrane" SubVolume1Ref="subdomain0" SubVolume2Ref="subdomain1" />
+        <SurfaceDescription NumSamplesX="101" NumSamplesY="101" NumSamplesZ="101" CutoffFrequency="0.3">
+          <VolumeRegion Name="subdomain00" RegionID="0" SubVolume="subdomain0" Size="732.257" Unit="um3" />
+          <VolumeRegion Name="subdomain11" RegionID="1" SubVolume="subdomain1" Size="267.743" Unit="um3" />
+          <MembraneRegion Name="membrane_subdomain00_subdomain11" VolumeRegion1="subdomain11" VolumeRegion2="subdomain00" Size="200.75261731312597" Unit="um2" />
+        </SurfaceDescription>
+      </Geometry>
+      <GeometryContext>
+        <FeatureMapping Feature="c0" GeometryClass="subdomain0" SubVolume="subdomain0" Size="50000.0" VolumePerUnitVolume="1.0">
+          <BoundariesTypes Xm="Flux" Xp="Flux" Ym="Flux" Yp="Flux" Zm="Flux" Zp="Flux" />
+        </FeatureMapping>
+        <FeatureMapping Feature="c1" GeometryClass="subdomain1" SubVolume="subdomain1" Size="50000.0" VolumePerUnitVolume="1.0">
+          <BoundariesTypes Xm="Flux" Xp="Flux" Ym="Flux" Yp="Flux" Zm="Flux" Zp="Flux" />
+        </FeatureMapping>
+        <MembraneMapping Membrane="m0" Size="6563.0" AreaPerUnitArea="1.0" CalculateVoltage="false" SpecificCapacitance="1.0" InitialVoltage="0.0" GeometryClass="subdomain0_subdomain1_membrane">
+          <BoundariesTypes Xm="Flux" Xp="Flux" Ym="Flux" Yp="Flux" Zm="Flux" Zp="Flux" />
+        </MembraneMapping>
+      </GeometryContext>
+      <ReactionContext>
+        <LocalizedCompoundSpec LocalizedCompoundRef="s0" ForceConstant="false" WellMixed="false" ForceContinuous="false">
+          <InitialConcentration>(1.0 + sin(x))</InitialConcentration>
+          <Diffusion>1.0E-4</Diffusion>
+          <Boundaries Xm="0.0" Xp="0.0" Ym="0.0" Yp="0.0" Zm="0.0" Zp="0.0" />
+        </LocalizedCompoundSpec>
+        <LocalizedCompoundSpec LocalizedCompoundRef="s1" ForceConstant="false" WellMixed="false" ForceContinuous="false">
+          <InitialConcentration>(1.0 + cos(x))</InitialConcentration>
+          <Diffusion>1.0E-4</Diffusion>
+          <Boundaries Xm="0.0" Xp="0.0" Ym="0.0" Yp="0.0" Zm="0.0" Zp="0.0" />
+        </LocalizedCompoundSpec>
+        <LocalizedCompoundSpec LocalizedCompoundRef="s3" ForceConstant="false" WellMixed="false" ForceContinuous="false">
+          <InitialConcentration>(1.0 + (sin(x) * cos(y)))</InitialConcentration>
+          <Diffusion>1.0E-4</Diffusion>
+          <Boundaries Xm="0.0" Xp="0.0" Ym="0.0" Yp="0.0" Zm="0.0" Zp="0.0" />
+        </LocalizedCompoundSpec>
+        <LocalizedCompoundSpec LocalizedCompoundRef="s2" ForceConstant="false" WellMixed="false" ForceContinuous="false">
+          <InitialConcentration>(1.0 + cos(y))</InitialConcentration>
+          <Diffusion>1.0000000000000002E-6</Diffusion>
+          <Boundaries Xm="0.0" Xp="0.0" Ym="0.0" Yp="0.0" Zm="0.0" Zp="0.0" />
+        </LocalizedCompoundSpec>
+        <ReactionSpec ReactionStepRef="r0" ReactionMapping="included" />
+        <ReactionSpec ReactionStepRef="r1" ReactionMapping="included" />
+        <ReactionSpec ReactionStepRef="r2" ReactionMapping="included" />
+      </ReactionContext>
+      <SpatialObjects>
+        <SpatialObject Name="vobj_subdomain00" Type="Volume" subVolume="subdomain0" regionId="0">
+          <QuantityCategoryList>
+            <QuantityCategory Name="VolumeCentroid" Enabled="false" />
+            <QuantityCategory Name="InteriorVelocity" Enabled="false" />
+            <QuantityCategory Name="VolumeRegionSize" Enabled="true" />
+          </QuantityCategoryList>
+        </SpatialObject>
+        <SpatialObject Name="vobj_subdomain11" Type="Volume" subVolume="subdomain1" regionId="1">
+          <QuantityCategoryList>
+            <QuantityCategory Name="VolumeCentroid" Enabled="false" />
+            <QuantityCategory Name="InteriorVelocity" Enabled="false" />
+            <QuantityCategory Name="VolumeRegionSize" Enabled="true" />
+          </QuantityCategoryList>
+        </SpatialObject>
+        <SpatialObject Name="sobj_subdomain11_subdomain00" Type="Surface" subVolumeInside="subdomain1" regionIdInside="1" subVolumeOutside="subdomain0" regionIdOutside="0">
+          <QuantityCategoryList>
+            <QuantityCategory Name="SurfaceNormal" Enabled="false" />
+            <QuantityCategory Name="SurfaceVelocity" Enabled="false" />
+            <QuantityCategory Name="DistanceToSurface" Enabled="false" />
+            <QuantityCategory Name="DirectionToSurface" Enabled="false" />
+            <QuantityCategory Name="SurfaceSize" Enabled="true" />
+          </QuantityCategoryList>
+        </SpatialObject>
+      </SpatialObjects>
+      <MicroscopeMeasurement Name="fluor">
+        <ConvolutionKernel Type="ProjectionZKernel" />
+      </MicroscopeMeasurement>
+    </SimulationSpec>
+    <pathwayModel>
+      <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:bp="http://www.biopax.org/release/biopax-level3.owl#" version="3.0" />
+    </pathwayModel>
+    <relationshipModel>
+      <RMNS version="3.0" />
+    </relationshipModel>
+    <vcmetadata>
+      <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" />
+      <nonrdfAnnotationList />
+      <uriBindingList />
+    </vcmetadata>
+  </BioModel>
+</vcml>
+
diff --git a/tests/fixtures/data/TinySpatialProject_Application0.vcml b/tests/fixtures/data/TinySpatialProject_Application0.vcml
index 9328bdb..0facf6a 100644
--- a/tests/fixtures/data/TinySpatialProject_Application0.vcml
+++ b/tests/fixtures/data/TinySpatialProject_Application0.vcml
@@ -1,6 +1,6 @@
 <?xml version="1.0" encoding="UTF-8"?>
-<!--This biomodel was generated in VCML Version Alpha_Version_7.7.0_build_15-->
-<vcml xmlns="http://sourceforge.net/projects/vcell/vcml" Version="Alpha_Version_7.7.0_build_15">
+<!--This biomodel was generated in VCML Version Alpha_Version_7.7.0_build_17-->
+<vcml xmlns="http://sourceforge.net/projects/vcell/vcml" Version="Alpha_Version_7.7.0_build_17">
   <BioModel Name="TinySpacialProject_Application0">
     <Model Name="unnamed">
       <ModelParameters>
@@ -13,125 +13,53 @@
       <Compound Name="s1">
         <Annotation>s1</Annotation>
       </Compound>
-      <Compound Name="s2" />
-      <Compound Name="s3" />
       <Feature Name="c0" />
-      <Feature Name="c1" />
-      <Membrane MembraneVoltage="Voltage_m0" Name="m0" />
       <LocalizedCompound Name="s0" SbmlName="s0" CompoundRef="s0" Structure="c0" OverrideName="true" />
       <LocalizedCompound Name="s1" SbmlName="s1" CompoundRef="s1" Structure="c0" OverrideName="true" />
-      <LocalizedCompound Name="s2" CompoundRef="s2" Structure="m0" OverrideName="true" />
-      <LocalizedCompound Name="s3" CompoundRef="s3" Structure="c1" OverrideName="true" />
       <SimpleReaction Structure="c0" Name="r0" Reversible="true" FluxOption="MolecularOnly" SbmlName="r0">
         <Reactant LocalizedCompoundRef="s0" Stoichiometry="1" />
         <Product LocalizedCompoundRef="s1" Stoichiometry="1" />
         <Kinetics KineticsType="GeneralKinetics">
-          <Parameter Name="J" Role="reaction rate" Unit="umol.l-1.s-1">((Kf_r0 * s0) - (1.0 * Kr_r0 * s1))</Parameter>
+          <Parameter Name="J" Role="reaction rate" Unit="umol.l-1.s-1">((Kf_r0 * s0) - (Kr_r0 * s1))</Parameter>
         </Kinetics>
       </SimpleReaction>
-      <SimpleReaction Structure="m0" Name="r1" Reversible="true" FluxOption="MolecularOnly">
-        <Reactant LocalizedCompoundRef="s0" Stoichiometry="1" />
-        <Product LocalizedCompoundRef="s2" Stoichiometry="1" />
-        <Kinetics KineticsType="MassAction">
-          <Parameter Name="J" Role="reaction rate" Unit="umol.dm-2.s-1">((Kf * s0) - (Kr * s2))</Parameter>
-          <Parameter Name="I" Role="inward current density" Unit="pA.dm-2">0.0</Parameter>
-          <Parameter Name="netValence" Role="net charge valence" Unit="1">1.0</Parameter>
-          <Parameter Name="Kf" Role="forward rate constant" Unit="l.dm-2.s-1">1.0</Parameter>
-          <Parameter Name="Kr" Role="reverse rate constant" Unit="s-1">1.0</Parameter>
-        </Kinetics>
-      </SimpleReaction>
-      <SimpleReaction Structure="m0" Name="r2" Reversible="true" FluxOption="MolecularOnly">
-        <Reactant LocalizedCompoundRef="s2" Stoichiometry="1" />
-        <Product LocalizedCompoundRef="s3" Stoichiometry="1" />
-        <Kinetics KineticsType="MassAction">
-          <Parameter Name="J" Role="reaction rate" Unit="umol.dm-2.s-1">((Kf * s2) - (Kr * s3))</Parameter>
-          <Parameter Name="I" Role="inward current density" Unit="pA.dm-2">0.0</Parameter>
-          <Parameter Name="netValence" Role="net charge valence" Unit="1">1.0</Parameter>
-          <Parameter Name="Kf" Role="forward rate constant" Unit="s-1">1.0</Parameter>
-          <Parameter Name="Kr" Role="reverse rate constant" Unit="l.dm-2.s-1">1.0</Parameter>
-        </Kinetics>
-      </SimpleReaction>
-      <Diagram Name="c1" Structure="c1">
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s0" LocationX="14" LocationY="34" />
-        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r0" LocationX="244" LocationY="128" />
-        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r2" LocationX="23" LocationY="88" />
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s3" LocationX="121" LocationY="84" />
-        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r1" LocationX="86" LocationY="49" />
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s2" LocationX="41" LocationY="81" />
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s1" LocationX="470" LocationY="220" />
-      </Diagram>
-      <Diagram Name="m0" Structure="m0">
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s0" LocationX="14" LocationY="34" />
-        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r0" LocationX="244" LocationY="128" />
-        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r2" LocationX="23" LocationY="88" />
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s3" LocationX="121" LocationY="84" />
-        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r1" LocationX="86" LocationY="49" />
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s2" LocationX="41" LocationY="81" />
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s1" LocationX="470" LocationY="220" />
-      </Diagram>
       <Diagram Name="c0" Structure="c0">
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s0" LocationX="14" LocationY="34" />
-        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r0" LocationX="244" LocationY="128" />
-        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r2" LocationX="23" LocationY="88" />
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s3" LocationX="121" LocationY="84" />
-        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r1" LocationX="86" LocationY="49" />
-        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s2" LocationX="41" LocationY="81" />
         <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s1" LocationX="470" LocationY="220" />
+        <SimpleReactionShape NodeReferenceModeAttrTag="full" SimpleReactionRef="r0" LocationX="244" LocationY="129" />
+        <LocalizedCompoundShape NodeReferenceModeAttrTag="full" LocalizedCompoundRef="s0" LocationX="14" LocationY="34" />
       </Diagram>
       <ModelUnitSystem VolumeSubstanceUnit="umol" MembraneSubstanceUnit="umol" LumpedReactionSubstanceUnit="umol" VolumeUnit="l" AreaUnit="dm2" LengthUnit="dm" TimeUnit="s" />
     </Model>
-    <SimulationSpec Name="unnamed_spatialGeom" Stochastic="false" UseConcentration="true" SpringSaLaD="false" RuleBased="false" MassConservationModelReduction="false" InsufficientIterations="false" InsufficientMaxMolecules="false" CharacteristicSize="0.10403160765672961">
+    <SimulationSpec Name="unnamed_spatialGeom" Stochastic="false" UseConcentration="true" SpringSaLaD="false" RuleBased="false" MassConservationModelReduction="false" InsufficientIterations="false" InsufficientMaxMolecules="false" CharacteristicSize="0.021090187684898395">
       <NetworkConstraints RbmMaxIteration="1" RbmMaxMoleculesPerSpecies="10" RbmSpeciesLimit="800" RbmReactionsLimit="2500" />
       <Annotation />
-      <Geometry Name="Geometry3" Dimension="3">
-        <Extent X="10.0" Y="10.0" Z="10.0" />
+      <Geometry Name="spatialGeom" Dimension="1">
+        <Extent X="10.0" Y="1.0" Z="1.0" />
         <Origin X="0.0" Y="0.0" Z="0.0" />
-        <SubVolume Name="subdomain1" Handle="1" Type="Analytical">
-          <AnalyticExpression>((((-5.0 + x) ^ 2.0) + ((-5.0 + y) ^ 2.0) + ((-5.0 + z) ^ 2.0)) &lt; 16.0)</AnalyticExpression>
-        </SubVolume>
         <SubVolume Name="subdomain0" Handle="0" Type="Analytical">
           <AnalyticExpression>1.0</AnalyticExpression>
         </SubVolume>
-        <SurfaceClass Name="subdomain0_subdomain1_membrane" SubVolume1Ref="subdomain0" SubVolume2Ref="subdomain1" />
-        <SurfaceDescription NumSamplesX="101" NumSamplesY="101" NumSamplesZ="101" CutoffFrequency="0.3">
-          <VolumeRegion Name="subdomain00" RegionID="0" SubVolume="subdomain0" Size="732.257" Unit="um3" />
-          <VolumeRegion Name="subdomain11" RegionID="1" SubVolume="subdomain1" Size="267.743" Unit="um3" />
-          <MembraneRegion Name="membrane_subdomain00_subdomain11" VolumeRegion1="subdomain11" VolumeRegion2="subdomain00" Size="200.75261731312597" Unit="um2" />
+        <SurfaceDescription NumSamplesX="50" NumSamplesY="1" NumSamplesZ="1" CutoffFrequency="0.3">
+          <VolumeRegion Name="subdomain00" RegionID="0" SubVolume="subdomain0" Size="10.0" Unit="um" />
         </SurfaceDescription>
       </Geometry>
       <GeometryContext>
         <FeatureMapping Feature="c0" GeometryClass="subdomain0" SubVolume="subdomain0" Size="50000.0" VolumePerUnitVolume="1.0">
           <BoundariesTypes Xm="Flux" Xp="Flux" Ym="Flux" Yp="Flux" Zm="Flux" Zp="Flux" />
         </FeatureMapping>
-        <FeatureMapping Feature="c1" GeometryClass="subdomain1" SubVolume="subdomain1" Size="50000.0" VolumePerUnitVolume="1.0">
-          <BoundariesTypes Xm="Flux" Xp="Flux" Ym="Flux" Yp="Flux" Zm="Flux" Zp="Flux" />
-        </FeatureMapping>
-        <MembraneMapping Membrane="m0" Size="6563.0" AreaPerUnitArea="1.0" CalculateVoltage="false" SpecificCapacitance="1.0" InitialVoltage="0.0" GeometryClass="subdomain0_subdomain1_membrane">
-          <BoundariesTypes Xm="Flux" Xp="Flux" Ym="Flux" Yp="Flux" Zm="Flux" Zp="Flux" />
-        </MembraneMapping>
       </GeometryContext>
       <ReactionContext>
         <LocalizedCompoundSpec LocalizedCompoundRef="s0" ForceConstant="false" WellMixed="false" ForceContinuous="false">
-          <InitialConcentration>(1.0 + sin(x))</InitialConcentration>
-          <Diffusion>1.0E-4</Diffusion>
+          <InitialConcentration>(100000.0 * x)</InitialConcentration>
+          <Diffusion>1.0E-9</Diffusion>
           <Boundaries Xm="0.0" Xp="0.0" />
         </LocalizedCompoundSpec>
         <LocalizedCompoundSpec LocalizedCompoundRef="s1" ForceConstant="false" WellMixed="false" ForceContinuous="false">
-          <InitialConcentration>(1.0 + cos(x))</InitialConcentration>
-          <Diffusion>1.0E-4</Diffusion>
+          <InitialConcentration>(10.0 - (1.0 * 100000.0 * x))</InitialConcentration>
+          <Diffusion>1.0E-9</Diffusion>
           <Boundaries Xm="0.0" Xp="0.0" />
         </LocalizedCompoundSpec>
-        <LocalizedCompoundSpec LocalizedCompoundRef="s3" ForceConstant="false" WellMixed="false" ForceContinuous="false">
-          <InitialConcentration>(1.0 + (sin(x) * cos(y)))</InitialConcentration>
-          <Diffusion>1.0E-4</Diffusion>
-        </LocalizedCompoundSpec>
-        <LocalizedCompoundSpec LocalizedCompoundRef="s2" ForceConstant="false" WellMixed="false" ForceContinuous="false">
-          <InitialConcentration>(1.0 + cos(y))</InitialConcentration>
-          <Diffusion>1.0000000000000002E-6</Diffusion>
-        </LocalizedCompoundSpec>
         <ReactionSpec ReactionStepRef="r0" ReactionMapping="included" />
-        <ReactionSpec ReactionStepRef="r1" ReactionMapping="included" />
-        <ReactionSpec ReactionStepRef="r2" ReactionMapping="included" />
       </ReactionContext>
       <MathDescription Name="unnamed_spatialGeom_generated">
         <Constant Name="_F_">96485.3321</Constant>
@@ -141,67 +69,31 @@
         <Constant Name="_PI_">3.141592653589793</Constant>
         <Constant Name="_R_">8314.46261815</Constant>
         <Constant Name="_T_">300.0</Constant>
-        <Constant Name="AreaPerUnitArea_m0">1.0</Constant>
         <Constant Name="K_millivolts_per_volt">1000.0</Constant>
         <Constant Name="Kf_r0">1.0</Constant>
-        <Constant Name="Kf_r1">1.0</Constant>
-        <Constant Name="Kf_r2">1.0</Constant>
         <Constant Name="KMOLE">0.001660538783162726</Constant>
         <Constant Name="Kr_r0">0.5</Constant>
-        <Constant Name="Kr_r1">1.0</Constant>
-        <Constant Name="Kr_r2">1.0</Constant>
-        <Constant Name="netValence_r1">1.0</Constant>
-        <Constant Name="netValence_r2">1.0</Constant>
         <Constant Name="s0_boundaryXm">0.0</Constant>
         <Constant Name="s0_boundaryXp">0.0</Constant>
-        <Constant Name="s0_diffusionRate">1.0E-4</Constant>
+        <Constant Name="s0_diffusionRate">1.0E-9</Constant>
         <Constant Name="s1_boundaryXm">0.0</Constant>
         <Constant Name="s1_boundaryXp">0.0</Constant>
-        <Constant Name="s1_diffusionRate">1.0E-4</Constant>
-        <Constant Name="s2_diffusionRate">1.0000000000000002E-6</Constant>
-        <Constant Name="s3_diffusionRate">1.0E-4</Constant>
-        <Constant Name="Voltage_m0">0.0</Constant>
+        <Constant Name="s1_diffusionRate">1.0E-9</Constant>
         <Constant Name="VolumePerUnitVolume_c0">1.0</Constant>
-        <Constant Name="VolumePerUnitVolume_c1">1.0</Constant>
         <VolumeVariable Name="s0" Domain="subdomain0" />
         <VolumeVariable Name="s1" Domain="subdomain0" />
-        <MembraneVariable Name="s2" Domain="subdomain0_subdomain1_membrane" />
-        <VolumeVariable Name="s3" Domain="subdomain1" />
-        <Function Name="J_r0" Domain="subdomain0">((Kf_r0 * s0) - (1.0 * Kr_r0 * s1))</Function>
-        <Function Name="J_r1" Domain="subdomain0_subdomain1_membrane">((Kf_r1 * s0) - (Kr_r1 * s2))</Function>
-        <Function Name="J_r2" Domain="subdomain0_subdomain1_membrane">((Kf_r2 * s2) - (Kr_r2 * s3))</Function>
-        <Function Name="KFlux_m0_c0" Domain="subdomain0_subdomain1_membrane">(AreaPerUnitArea_m0 / VolumePerUnitVolume_c0)</Function>
-        <Function Name="KFlux_m0_c1" Domain="subdomain0_subdomain1_membrane">(AreaPerUnitArea_m0 / VolumePerUnitVolume_c1)</Function>
-        <Function Name="s0_init_umol_l_1" Domain="subdomain0">(1.0 + sin(x))</Function>
-        <Function Name="s1_init_umol_l_1" Domain="subdomain0">(1.0 + cos(x))</Function>
-        <Function Name="s2_init_umol_dm_2" Domain="subdomain0_subdomain1_membrane">(1.0 + cos(y))</Function>
-        <Function Name="s3_init_umol_l_1" Domain="subdomain1">(1.0 + (sin(x) * cos(y)))</Function>
+        <Function Name="J_r0" Domain="subdomain0">((Kf_r0 * s0) - (Kr_r0 * s1))</Function>
+        <Function Name="s0_init_umol_l_1" Domain="subdomain0">(100000.0 * x)</Function>
+        <Function Name="s1_init_umol_l_1" Domain="subdomain0">(10.0 - (1.0 * 100000.0 * x))</Function>
         <Function Name="Size_c0" Domain="subdomain0">(VolumePerUnitVolume_c0 * vcRegionVolume('subdomain0'))</Function>
-        <Function Name="Size_c1" Domain="subdomain1">(VolumePerUnitVolume_c1 * vcRegionVolume('subdomain1'))</Function>
-        <Function Name="Size_m0" Domain="subdomain0_subdomain1_membrane">(AreaPerUnitArea_m0 * vcRegionArea('subdomain0_subdomain1_membrane'))</Function>
-        <Function Name="sobj_subdomain11_subdomain00_size" Domain="subdomain0_subdomain1_membrane">vcRegionArea('subdomain0_subdomain1_membrane')</Function>
         <Function Name="vobj_subdomain00_size" Domain="subdomain0">vcRegionVolume('subdomain0')</Function>
-        <Function Name="vobj_subdomain11_size" Domain="subdomain1">vcRegionVolume('subdomain1')</Function>
-        <CompartmentSubDomain Name="subdomain1">
-          <BoundaryType Boundary="Xm" Type="Flux" />
-          <BoundaryType Boundary="Xp" Type="Flux" />
-          <BoundaryType Boundary="Ym" Type="Flux" />
-          <BoundaryType Boundary="Yp" Type="Flux" />
-          <BoundaryType Boundary="Zm" Type="Flux" />
-          <BoundaryType Boundary="Zp" Type="Flux" />
-          <PdeEquation Name="s3" SolutionType="Unknown">
-            <Rate>0.0</Rate>
-            <Diffusion>s3_diffusionRate</Diffusion>
-            <Initial>s3_init_umol_l_1</Initial>
-          </PdeEquation>
-        </CompartmentSubDomain>
         <CompartmentSubDomain Name="subdomain0">
           <BoundaryType Boundary="Xm" Type="Flux" />
           <BoundaryType Boundary="Xp" Type="Flux" />
-          <BoundaryType Boundary="Ym" Type="Flux" />
-          <BoundaryType Boundary="Yp" Type="Flux" />
-          <BoundaryType Boundary="Zm" Type="Flux" />
-          <BoundaryType Boundary="Zp" Type="Flux" />
+          <BoundaryType Boundary="Ym" Type="Value" />
+          <BoundaryType Boundary="Yp" Type="Value" />
+          <BoundaryType Boundary="Zm" Type="Value" />
+          <BoundaryType Boundary="Zp" Type="Value" />
           <PdeEquation Name="s0" SolutionType="Unknown">
             <Boundaries Xm="s0_boundaryXm" Xp="s0_boundaryXp" />
             <Rate>- J_r0</Rate>
@@ -215,38 +107,13 @@
             <Initial>s1_init_umol_l_1</Initial>
           </PdeEquation>
         </CompartmentSubDomain>
-        <MembraneSubDomain Name="subdomain0_subdomain1_membrane" InsideCompartment="subdomain0" OutsideCompartment="subdomain1">
-          <BoundaryType Boundary="Xm" Type="Value" />
-          <BoundaryType Boundary="Xp" Type="Value" />
-          <BoundaryType Boundary="Ym" Type="Value" />
-          <BoundaryType Boundary="Yp" Type="Value" />
-          <BoundaryType Boundary="Zm" Type="Value" />
-          <BoundaryType Boundary="Zp" Type="Value" />
-          <PdeEquation Name="s2" SolutionType="Unknown">
-            <Rate>(J_r1 - J_r2)</Rate>
-            <Diffusion>s2_diffusionRate</Diffusion>
-            <Initial>s2_init_umol_dm_2</Initial>
-          </PdeEquation>
-          <JumpCondition Name="s0">
-            <InFlux>( - 1.0 * KFlux_m0_c0 * J_r1)</InFlux>
-            <OutFlux>0.0</OutFlux>
-          </JumpCondition>
-          <JumpCondition Name="s1">
-            <InFlux>0.0</InFlux>
-            <OutFlux>0.0</OutFlux>
-          </JumpCondition>
-          <JumpCondition Name="s3">
-            <InFlux>0.0</InFlux>
-            <OutFlux>(KFlux_m0_c1 * J_r2)</OutFlux>
-          </JumpCondition>
-        </MembraneSubDomain>
       </MathDescription>
       <Simulation Name="Simulation0">
         <SolverTaskDescription TaskType="Unsteady" UseSymbolicJacobian="false" Solver="Sundials Stiff PDE Solver (Variable Time Step)">
-          <TimeBound StartTime="0.0" EndTime="5.0" />
+          <TimeBound StartTime="0.0" EndTime="1.0" />
           <TimeStep DefaultTime="0.05" MinTime="0.0" MaxTime="0.1" />
           <ErrorTolerance Absolut="1.0E-9" Relative="1.0E-7" />
-          <OutputOptions OutputTimeStep="0.1" />
+          <OutputOptions OutputTimeStep="0.05" />
           <SundialsSolverOptions>
             <maxOrderAdvection>2</maxOrderAdvection>
           </SundialsSolverOptions>
@@ -254,7 +121,7 @@
         </SolverTaskDescription>
         <MathOverrides />
         <MeshSpecification>
-          <Size X="51" Y="51" Z="51" />
+          <Size X="200" Y="1" Z="1" />
         </MeshSpecification>
       </Simulation>
       <SpatialObjects>
@@ -265,22 +132,6 @@
             <QuantityCategory Name="VolumeRegionSize" Enabled="true" />
           </QuantityCategoryList>
         </SpatialObject>
-        <SpatialObject Name="vobj_subdomain11" Type="Volume" subVolume="subdomain1" regionId="1">
-          <QuantityCategoryList>
-            <QuantityCategory Name="VolumeCentroid" Enabled="false" />
-            <QuantityCategory Name="InteriorVelocity" Enabled="false" />
-            <QuantityCategory Name="VolumeRegionSize" Enabled="true" />
-          </QuantityCategoryList>
-        </SpatialObject>
-        <SpatialObject Name="sobj_subdomain11_subdomain00" Type="Surface" subVolumeInside="subdomain1" regionIdInside="1" subVolumeOutside="subdomain0" regionIdOutside="0">
-          <QuantityCategoryList>
-            <QuantityCategory Name="SurfaceNormal" Enabled="false" />
-            <QuantityCategory Name="SurfaceVelocity" Enabled="false" />
-            <QuantityCategory Name="DistanceToSurface" Enabled="false" />
-            <QuantityCategory Name="DirectionToSurface" Enabled="false" />
-            <QuantityCategory Name="SurfaceSize" Enabled="true" />
-          </QuantityCategoryList>
-        </SpatialObject>
       </SpatialObjects>
       <MicroscopeMeasurement Name="fluor">
         <ConvolutionKernel Type="ProjectionZKernel" />
diff --git a/tests/fixtures/data/TinySpatialProject_Application0_from_python.vcml b/tests/fixtures/data/TinySpatialProject_Application0_from_python.vcml
new file mode 100644
index 0000000..bd032c1
--- /dev/null
+++ b/tests/fixtures/data/TinySpatialProject_Application0_from_python.vcml
@@ -0,0 +1,70 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<vcml Version="Alpha_Version_7.7.0_build_15">
+  <BioModel Name="TinySpacialProject_Application0">
+    <Model Name="unnamed">
+      <ModelParameters>
+        <Parameter Name="Kf_r0" Role="user defined" Unit="s-1">1.0</Parameter>
+        <Parameter Name="Kr_r0" Role="user defined" Unit="s-1">0.5</Parameter>
+      </ModelParameters>
+      <Compound Name="s0">
+        <Annotation>s0</Annotation>
+      </Compound>
+      <Compound Name="s1">
+        <Annotation>s1</Annotation>
+      </Compound>
+      <Feature Name="c0"/>
+      <LocalizedCompound Name="s0" CompoundRef="s0" Structure="c0"/>
+      <LocalizedCompound Name="s1" CompoundRef="s1" Structure="c0"/>
+      <SimpleReaction Structure="c0" Name="r0">
+        <Reactant LocalizedCompoundRef="s0" Stoichiometry="1"/>
+        <Product LocalizedCompoundRef="s1" Stoichiometry="1"/>
+        <Kinetics KineticsType="GeneralKinetics">
+          <Parameter Name="J" Role="reaction rate" Unit="umol.l-1.s-1">((Kf_r0 * s0) - (Kr_r0 * s1))</Parameter>
+        </Kinetics>
+      </SimpleReaction>
+    </Model>
+    <SimulationSpec Name="unnamed_spatialGeom">
+      <Geometry Name="spatialGeom" Dimension="1">
+        <Extent X="10.0" Y="1.0" Z="1.0"/>
+        <Origin X="0.0" Y="0.0" Z="0.0"/>
+        <SubVolume Name="subdomain0" Handle="0" Type="Analytical">
+          <AnalyticExpression>1.0</AnalyticExpression>
+        </SubVolume>
+      </Geometry>
+      <GeometryContext>
+        <FeatureMapping Feature="c0" GeometryClass="subdomain0" VolumePerUnitVolume="1.0">
+          <BoundariesTypes Xm="Flux" Xp="Flux" Ym="Flux" Yp="Flux" Zm="Flux" Zp="Flux"/>
+        </FeatureMapping>
+      </GeometryContext>
+      <ReactionContext>
+        <LocalizedCompoundSpec LocalizedCompoundRef="s0">
+          <InitialConcentration>(100000.0 * x)</InitialConcentration>
+          <Diffusion>1e-09</Diffusion>
+          <Boundaries Xm="0.0" Xp="0.0"/>
+        </LocalizedCompoundSpec>
+        <LocalizedCompoundSpec LocalizedCompoundRef="s1">
+          <InitialConcentration>(10.0 - (1.0 * 100000.0 * x))</InitialConcentration>
+          <Diffusion>1e-09</Diffusion>
+          <Boundaries Xm="0.0" Xp="0.0"/>
+        </LocalizedCompoundSpec>
+        <ReactionSpec ReactionStepRef="r0" ReactionMapping="included"/>
+      </ReactionContext>
+      <Simulation Name="Simulation0">
+        <SolverTaskDescription TaskType="Unsteady" UseSymbolicJacobian="false" Solver="Sundials Stiff PDE Solver (Variable Time Step)">
+          <TimeBound StartTime="0.0" EndTime="1.0"/>
+          <TimeStep DefaultTime="0.05" MinTime="0.0" MaxTime="0.1"/>
+          <ErrorTolerance Absolut="1.0E-9" Relative="1.0E-7"/>
+          <OutputOptions OutputTimeStep="0.05"/>
+          <SundialsSolverOptions>
+            <maxOrderAdvection>2</maxOrderAdvection>
+          </SundialsSolverOptions>
+          <NumberProcessors>1</NumberProcessors>
+        </SolverTaskDescription>
+        <MathOverrides/>
+        <MeshSpecification>
+          <Size X="200" Y="1" Z="1"/>
+        </MeshSpecification>
+      </Simulation>
+    </SimulationSpec>
+  </BioModel>
+</vcml>
diff --git a/tests/fixtures/vcell_model_fixtures.py b/tests/fixtures/vcell_model_fixtures.py
index 9f035f6..751a2d7 100644
--- a/tests/fixtures/vcell_model_fixtures.py
+++ b/tests/fixtures/vcell_model_fixtures.py
@@ -9,3 +9,8 @@
 @pytest.fixture
 def vcml_spatial_model_1d_path() -> Path:
     return FIXTURE_DATA_DIR / "TinySpatialProject_Application0.vcml"
+
+
+@pytest.fixture
+def vcml_spatial_small_3d_path() -> Path:
+    return FIXTURE_DATA_DIR / "SmallSpacialProject_3D.vcml"
diff --git a/tests/vcml/test_vcml_reader.py b/tests/vcml/test_vcml_reader.py
index 44d1b74..82ed4af 100644
--- a/tests/vcml/test_vcml_reader.py
+++ b/tests/vcml/test_vcml_reader.py
@@ -3,7 +3,7 @@
 import pyvcell.vcml as vc
 
 
-def test_vcml_reader(vcml_spatial_model_1d_path: Path) -> None:
+def test_vcml_reader_1D(vcml_spatial_model_1d_path: Path) -> None:
     assert vcml_spatial_model_1d_path.is_file()
 
     with open(vcml_spatial_model_1d_path) as f:
@@ -15,6 +15,61 @@ def test_vcml_reader(vcml_spatial_model_1d_path: Path) -> None:
     assert model is not None and model.name == "unnamed"
 
     assert [p.name for p in model.model_parameters] == ["Kf_r0", "Kr_r0"]
+    assert [r.name for r in model.reactions] == ["r0"]
+    assert [(c.name, c.dim) for c in model.compartments] == [("c0", 3)]
+    assert [(s.name, s.structure_name) for s in model.species] == [
+        ("s0", "c0"),
+        ("s1", "c0"),
+    ]
+
+    r0: vc.Reaction = model.reactions[0]
+    assert [(r.name, r.stoichiometry) for r in r0.reactants] == [("s0", 1)]
+    assert [(p.name, p.stoichiometry) for p in r0.products] == [("s1", 1)]
+    assert r0.kinetics is not None and r0.kinetics.kinetics_type == "GeneralKinetics"
+    assert r0.compartment_name == "c0"
+    assert {p.name: p.value for p in r0.kinetics.kinetics_parameters} == {"J": "((Kf_r0 * s0) - (Kr_r0 * s1))"}
+
+    assert [a.name for a in biomodel.applications] == ["unnamed_spatialGeom"]
+    app0 = biomodel.applications[0]
+    assert app0.stochastic is False
+
+    geom = app0.geometry
+    assert (geom.name, geom.dim) == ("spatialGeom", 1)
+    assert geom.extent == (10.0, 1.0, 1.0)
+    assert geom.origin == (0.0, 0.0, 0.0)
+    assert [(sv.name, sv.handle, sv.subvolume_type.name, sv.analytic_expr) for sv in geom.subvolumes] == [
+        ("subdomain0", 0, "analytic", "1.0")
+    ]
+    assert [(sc.name, sc.subvolume_ref_1, sc.subvolume_ref_2) for sc in geom.surface_classes] == []
+
+    assert [
+        (cm.compartment_name, cm.geometry_class_name, cm.unit_size, cm.boundary_types)
+        for cm in app0.compartment_mappings
+    ] == [("c0", "subdomain0", 1.0, ["flux", "flux", "flux", "flux", "flux", "flux"])]
+
+    assert [
+        (sm.species_name, sm.diffusion_coefficient, sm.initial_concentration, sm.boundary_values)
+        for sm in app0.species_mappings
+    ] == [
+        ("s0", 1e-09, "(100000.0 * x)", [0.0, 0.0, None, None, None, None]),
+        ("s1", 1e-09, "(10.0 - (1.0 * 100000.0 * x))", [0.0, 0.0, None, None, None, None]),
+    ]
+
+    assert [[(rm.reaction_name, rm.included) for rm in app0.reaction_mappings]] == [[("r0", True)]]
+
+
+def test_vcml_reader_3D(vcml_spatial_small_3d_path: Path) -> None:
+    assert vcml_spatial_small_3d_path.is_file()
+
+    with open(vcml_spatial_small_3d_path) as f:
+        xml_string = f.read()
+
+    biomodel = vc.VcmlReader.parse_biomodel(xml_string)
+    assert biomodel is not None and biomodel.name == "TinySpacialProject_Application0_unnamed_spatialGeom"
+    model = biomodel.model
+    assert model is not None and model.name == "unnamed"
+
+    assert [p.name for p in model.model_parameters] == ["Kf_r0", "Kr_r0", "Kf_r1", "Kr_r1", "Kf_r2", "Kr_r2"]
     assert [r.name for r in model.reactions] == ["r0", "r1", "r2"]
     assert [(c.name, c.dim) for c in model.compartments] == [("c0", 3), ("c1", 3), ("m0", 2)]
     assert [(s.name, s.structure_name) for s in model.species] == [
@@ -34,27 +89,23 @@ def test_vcml_reader(vcml_spatial_model_1d_path: Path) -> None:
     r1: vc.Reaction = model.reactions[1]
     assert [(r.name, r.stoichiometry) for r in r1.reactants] == [("s0", 1)]
     assert [(p.name, p.stoichiometry) for p in r1.products] == [("s2", 1)]
-    assert r1.kinetics is not None and r1.kinetics.kinetics_type == "MassAction"
+    assert r1.kinetics is not None and r1.kinetics.kinetics_type == "GeneralKinetics"
     assert r1.compartment_name == "m0"
     assert {p.name: p.value for p in r1.kinetics.kinetics_parameters} == {
-        "J": "((Kf * s0) - (Kr * s2))",
-        "I": 0,
-        "netValence": 1,
-        "Kf": 1,
-        "Kr": 1,
+        "I": 0.0,
+        "J": "((Kf_r1 * s0) - (Kr_r1 * s2))",
+        "netValence": 1.0,
     }
 
     r2: vc.Reaction = model.reactions[2]
     assert [(r.name, r.stoichiometry) for r in r2.reactants] == [("s2", 1)]
     assert [(p.name, p.stoichiometry) for p in r2.products] == [("s3", 1)]
-    assert r2.kinetics is not None and r2.kinetics.kinetics_type == "MassAction"
+    assert r2.kinetics is not None and r2.kinetics.kinetics_type == "GeneralKinetics"
     assert r2.compartment_name == "m0"
     assert {p.name: p.value for p in r2.kinetics.kinetics_parameters} == {
-        "J": "((Kf * s2) - (Kr * s3))",
-        "I": 0,
-        "netValence": 1,
-        "Kf": 1,
-        "Kr": 1,
+        "I": 0.0,
+        "J": "((Kf_r2 * s2) - (Kr_r2 * s3))",
+        "netValence": 1.0,
     }
 
     assert [a.name for a in biomodel.applications] == ["unnamed_spatialGeom"]
@@ -62,11 +113,11 @@ def test_vcml_reader(vcml_spatial_model_1d_path: Path) -> None:
     assert app0.stochastic is False
 
     geom = app0.geometry
-    assert (geom.name, geom.dim) == ("Geometry3", 3)
+    assert (geom.name, geom.dim) == ("spatialGeom", 3)
     assert geom.extent == (10.0, 10.0, 10.0)
     assert geom.origin == (0.0, 0.0, 0.0)
     assert [(sv.name, sv.handle, sv.subvolume_type.name, sv.analytic_expr) for sv in geom.subvolumes] == [
-        ("subdomain1", 1, "analytic", "((((-5.0 + x) ^ 2.0) + ((-5.0 + y) ^ 2.0) + ((-5.0 + z) ^ 2.0)) < 16.0)"),
+        ("subdomain1", 1, "analytic", "((pow((-5.0 + x),2.0) + pow((-5.0 + y),2.0) + pow((-5.0 + z),2.0)) < 16.0)"),
         ("subdomain0", 0, "analytic", "1.0"),
     ]
     assert [(sc.name, sc.subvolume_ref_1, sc.subvolume_ref_2) for sc in geom.surface_classes] == [
@@ -86,10 +137,10 @@ def test_vcml_reader(vcml_spatial_model_1d_path: Path) -> None:
         (sm.species_name, sm.diffusion_coefficient, sm.initial_concentration, sm.boundary_values)
         for sm in app0.species_mappings
     ] == [
-        ("s0", 0.0001, "(1.0 + sin(x))", [0.0, 0.0, None, None, None, None]),
-        ("s1", 0.0001, "(1.0 + cos(x))", [0.0, 0.0, None, None, None, None]),
-        ("s3", 0.0001, "(1.0 + (sin(x) * cos(y)))", []),
-        ("s2", 1.0000000000000002e-06, "(1.0 + cos(y))", []),
+        ("s0", 0.0001, "(1.0 + sin(x))", [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]),
+        ("s1", 0.0001, "(1.0 + cos(x))", [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]),
+        ("s3", 0.0001, "(1.0 + (sin(x) * cos(y)))", [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]),
+        ("s2", 1.0000000000000002e-06, "(1.0 + cos(y))", [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]),
     ]
 
     assert [[(rm.reaction_name, rm.included) for rm in app0.reaction_mappings]] == [
diff --git a/tests/vcml/test_vcml_writer.py b/tests/vcml/test_vcml_writer.py
index ee8db84..68f72eb 100644
--- a/tests/vcml/test_vcml_writer.py
+++ b/tests/vcml/test_vcml_writer.py
@@ -3,7 +3,7 @@
 import pyvcell.vcml as vc
 
 
-def test_vcml_writer(vcml_spatial_model_1d_path: Path) -> None:
+def test_vcml_writer_1D(vcml_spatial_model_1d_path: Path) -> None:
     assert vcml_spatial_model_1d_path.is_file()
 
     with open(vcml_spatial_model_1d_path) as f:
@@ -41,3 +41,43 @@ def test_vcml_writer(vcml_spatial_model_1d_path: Path) -> None:
         f.write(new_vcml_str)
 
     assert biomodel == new_biomodel
+
+
+def test_vcml_writer_3D(vcml_spatial_small_3d_path: Path) -> None:
+    assert vcml_spatial_small_3d_path.is_file()
+
+    with open(vcml_spatial_small_3d_path) as f:
+        orig_vcml_string = f.read()
+
+    biomodel = vc.VcmlReader.parse_biomodel(orig_vcml_string)
+    document = vc.VCMLDocument(biomodel=biomodel)
+
+    new_vcml_str: str = vc.VcmlWriter().write_vcml(document=document)
+    assert new_vcml_str is not None
+
+    new_biomodel = vc.VcmlReader.parse_biomodel(new_vcml_str)
+
+    assert biomodel is not None and new_biomodel is not None
+    assert biomodel.name == new_biomodel.name
+    assert biomodel.model is not None and new_biomodel.model is not None
+    assert biomodel.model.name == new_biomodel.model.name
+    assert biomodel.model.model_parameters == new_biomodel.model.model_parameters
+
+    for i, reaction in enumerate(biomodel.model.reactions):
+        new_reaction = new_biomodel.model.reactions[i]
+        assert reaction.name == new_reaction.name
+        assert reaction.reactants == new_reaction.reactants
+        assert reaction.products == new_reaction.products
+        assert reaction.kinetics == new_reaction.kinetics
+
+    for i, compartment in enumerate(biomodel.model.compartments):
+        new_compartment = new_biomodel.model.compartments[i]
+        assert compartment.name == new_compartment.name
+        assert compartment.dim == new_compartment.dim
+
+    # write out the new vcml string to a file so that it can be compared with original vcml file
+    new_vcml_path = vcml_spatial_small_3d_path.with_name("new_vcml_3d.xml")
+    with open(new_vcml_path, "w") as f:
+        f.write(new_vcml_str)
+
+    assert biomodel == new_biomodel