diff --git a/README.md b/README.md
index c58a9db..73eab01 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# Dynamic PET <img src="docs/logo.png" align="right" height="150" />
+# Dynamic PET <img src="docs/logo.png" align="right" height="150" alt="Dynamic PET logo"/>
 
 <!-- [![PyPI](https://img.shields.io/pypi/v/dynamicpet.svg)][pypi_]
 [![Status](https://img.shields.io/pypi/status/dynamicpet.svg)][status]
@@ -34,6 +34,7 @@ Methods implemented in the CLI include:
 
 - Denoising
   - HighlY constrained backPRojection method constraining the backprojections to Local Regions of interest ([HYPR-LR])
+  - Nonlocal EStimation of multispectral MAgnitudes ([NESMA])
 - Reference tissue-based modeling
   - Standardized Uptake Value Ratio (SUVR)
   - Logan Reference Tissue Model ([LRTM])
@@ -47,6 +48,7 @@ Several implementations of estimating SRTM parameters are available:
 
 [lrtm]: https://doi.org/10.1097/00004647-199609000-00008
 [hypr-lr]: https://doi.org/10.2967/jnumed.109.073999
+[nesma]: https://doi.org/10.1111/jon.12537
 
 ## Requirements
 
@@ -60,8 +62,8 @@ _Dynamic PET_ requires Python 3.11+ and the following modules:
 You can install _Dynamic PET_ via [pip] after cloning the repository:
 
 ```console
-$ git clone https://github.com/bilgelm/dynamicpet.git
-$ pip install -e dynamicpet
+git clone https://github.com/bilgelm/dynamicpet.git
+pip install -e dynamicpet
 ```
 
 ## Usage
@@ -78,13 +80,13 @@ You will then need to create a binary mask that is in the same space as the PET
 After installing _Dynamic PET_ as described above, execute:
 
 ```console
-$ kineticmodel PET --model SRTMZhou2003 --refmask <REFMASK> --outputdir <OUTPUTDIR> --fwhm 5
+kineticmodel PET --model SRTMZhou2003 --refmask <REFMASK> --outputdir <OUTPUTDIR> --fwhm 5
 ```
 
 where
 
 ```console
-$ PET=<OPENNEURODATA>/ds001705-download/sub-000101/ses-baseline/pet/sub-000101_ses-baseline_pet.nii
+PET=<OPENNEURODATA>/ds001705-download/sub-000101/ses-baseline/pet/sub-000101_ses-baseline_pet.nii
 ```
 
 Before running these commands, replace
diff --git a/pyproject.toml b/pyproject.toml
index 45e1dee..f14ef51 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "dynamicpet"
-version = "0.1.3"
+version = "0.1.4"
 description = "Dynamic PET"
 authors = ["Murat Bilgel <bilgelm@gmail.com>"]
 license = "MIT"
diff --git a/src/dynamicpet/__init__.py b/src/dynamicpet/__init__.py
index 00d86de..7e694ff 100644
--- a/src/dynamicpet/__init__.py
+++ b/src/dynamicpet/__init__.py
@@ -1 +1,6 @@
 """Dynamic PET."""
+
+import importlib.metadata
+
+
+__version__ = importlib.metadata.version("dynamicpet")
diff --git a/src/dynamicpet/__main__.py b/src/dynamicpet/__main__.py
index daab838..61b286e 100644
--- a/src/dynamicpet/__main__.py
+++ b/src/dynamicpet/__main__.py
@@ -2,7 +2,9 @@
 
 import csv
 import os
+import re
 import warnings
+from json import dump as json_dump
 
 import click
 import numpy as np
@@ -10,6 +12,7 @@
 from nibabel.loadsave import load as nib_load
 from nibabel.spatialimages import SpatialImage
 
+from dynamicpet import __version__
 from dynamicpet.denoise import hypr
 from dynamicpet.denoise import nesma
 from dynamicpet.kineticmodel.kineticmodel import KineticModel
@@ -100,7 +103,7 @@ def denoise(
 ) -> None:
     """Perform dynamic PET denoising.
 
-    Outputs will have a '_<method>' suffix.
+    Outputs will have a '_desc-<method>_pet' suffix.
 
     PET: 4-D PET image
     """
@@ -108,11 +111,23 @@ def denoise(
     pet_img = petbidsimage_load(pet, json)
 
     if method == "HYPRLR":
+        doc = hypr.hypr_lr.__doc__
+
+        # parameters not relevant to HYPR-LR
+        mask = None
+        window_half_size = None
+        thresh = None
+
         if fwhm:
             res = hypr.hypr_lr(pet_img, fwhm)
         else:
             raise ValueError("fwhm must be specified for HYPR-LR")
     elif method == "NESMA":
+        doc = nesma.nesma_semiadaptive.__doc__
+
+        # parameter not relevant to NESMA
+        fwhm = None
+
         if mask:
             mask_img: SpatialImage = nib_load(mask)  # type: ignore
             # check that mask is in the same space as pet
@@ -138,10 +153,37 @@ def denoise(
     if outputdir:
         os.makedirs(outputdir, exist_ok=True)
         bname = os.path.basename(froot)
-        froot = os.path.join(outputdir, bname)
-    output = froot + "_" + method.lower() + ext + addext
+        # if the input file name follows the PET-BIDS convention, it should end
+        # with "_pet". Need to move this to the end of the new file name to
+        # maintain compatibility with the PET-BIDS Derivatives convention.
+        froot = re.sub("_pet$", "", os.path.join(outputdir, bname))
+    output = froot + "_desc-" + method.lower() + "_pet" + ext + addext
+    output_json = froot + "_desc-" + method.lower() + "_pet.json"
     res.to_filename(output)
 
+    cmd = (
+        f"denoise --method {method} "
+        + (f"--fwhm {fwhm} " if fwhm else "")
+        + (f"--mask {mask} " if mask else "")
+        + (f"--window_half_size {window_half_size} " if window_half_size else "")
+        + (f"--thresh {thresh} " if thresh else "")
+        + (f"--outputdir {outputdir} " if outputdir else "")
+        + (f"--json {json} " if json else "")
+        + pet
+    )
+    derivative_json_dict = {
+        "Description": (
+            re.sub(r"\s+", " ", doc.split("Args:")[0]).strip() if doc else ""
+        ),
+        # "Sources": [pet],
+        "SoftwareName": "dynamicpet",
+        "SoftwareVersion": __version__,
+        "CommandLine": cmd,
+    }
+
+    with open(output_json, "w") as f:
+        json_dump(derivative_json_dict, f, indent=4)
+
 
 @click.command()
 @click.argument("pet", type=str)
@@ -200,9 +242,11 @@ def denoise(
     ),
 )
 @click.option(
-    "--start", default=None, type=float, help="Start of time window for model"
+    "--start", default=None, type=float, help="Start of time window for model in min"
+)
+@click.option(
+    "--end", default=None, type=float, help="End of time window for model in min"
 )
-@click.option("--end", default=None, type=float, help="End of time window for model")
 @click.option(
     "--fwhm",
     default=None,
@@ -226,7 +270,7 @@ def denoise(
         "'rect' is rectangular integration."
     ),
 )
-def kineticmodel(
+def kineticmodel(  # noqa: C901
     pet: str,
     model: str,
     refroi: str | None,
@@ -242,7 +286,7 @@ def kineticmodel(
 ) -> None:
     """Fit a reference tissue model to a dynamic PET image or TACs.
 
-    Outputs will have a '_km-<model>_kp-<parameter>' suffix.
+    Outputs will have a '_model-<model>_meas-<parameter>' suffix.
 
     PET: 4-D PET image (can be 3-D if model is SUVR) or a 2-D tabular TACs tsv file
     """
@@ -262,16 +306,27 @@ def kineticmodel(
         pet_img = pet_img.extract(start, end)
         reftac = reftac.extract(start, end)
 
+    if fwhm and model not in ["srtmzhou2003"]:
+        fwhm = None
+        warnings.warn(
+            "--fwhm argument is not relevant for this model, will be ignored",
+            RuntimeWarning,
+            stacklevel=2,
+        )
+
     # fit kinetic model
     km: KineticModel
     match model:
         case "suvr":
+            model_abbr = "SUVR"
             km = SUVR(reftac, pet_img)
             km.fit(mask=petmask_img_mat)
         case "srtmlammertsma1996":
+            model_abbr = "SRTM"
             km = SRTMLammertsma1996(reftac, pet_img)
             km.fit(mask=petmask_img_mat, weight_by=weight_by)
         case "srtmzhou2003":
+            model_abbr = "SRTM"
             km = SRTMZhou2003(reftac, pet_img)
             km.fit(
                 mask=petmask_img_mat,
@@ -288,7 +343,13 @@ def kineticmodel(
         froot = os.path.join(outputdir, bname)
 
     if isinstance(pet_img, PETBIDSMatrix):
-        output = froot + "_km-" + model.replace(".", "") + ext
+        # if the input file name follows the PET-BIDS Derivatives convention,
+        # it should end with "_tacs". Need to remove this to maintain
+        # compatibility with the PET-BIDS Derivatives convention.
+        froot = re.sub("_tacs$", "", froot)
+
+        output = froot + "_model-" + model_abbr + "_kinpar" + ext
+        output_json = froot + "_model-" + model_abbr + "_kinpar.json"
         data = np.empty((len(km.parameters), pet_img.num_elements))
         for i, param in enumerate(km.parameters.keys()):
             data[i] = km.get_parameter(param)
@@ -299,16 +360,73 @@ def kineticmodel(
             for i, elem in enumerate(pet_img.elem_names):
                 tsvwriter.writerow([elem] + datat[i].tolist())
     else:
+        # if the input file name follows the PET-BIDS convention, it should end
+        # with "_pet". Need to remove this to maintain compatibility with the
+        # PET-BIDS Derivatives convention.
+        froot = re.sub("_pet$", "", froot)
+
         # save estimated parameters as image
         for param in km.parameters.keys():
             res_img: SpatialImage = km.get_parameter(param)  # type: ignore
             output = (
-                froot + "_km-" + model.replace(".", "") + "_kp-" + param + ext + addext
+                froot
+                + "_model-"
+                + model_abbr
+                + "_meas-"
+                + param
+                + "_mimap"
+                + ext
+                + addext
             )
             res_img.to_filename(output)
-
-    # also need to save a json PET BIDS derivative file
-    # TODO
+        output_json = froot + "_model-" + model_abbr + "_mimap.json"
+
+    # save json PET BIDS derivative file
+    inputvalues = [start, end]
+    inputvalueslabels = [
+        "Start of time window for model",
+        "End of time window for model",
+    ]
+    inputvaluesunits = ["min", "min"]
+
+    if fwhm:
+        inputvalues += [fwhm]
+        inputvalueslabels += ["Full width at half max"]
+        inputvaluesunits += ["mm"]
+
+    cmd = (
+        f"kineticmodel --model {model} "
+        + (f"--refroi {refroi} " if refroi else f"--refmask {refmask} ")
+        + (f"--outputdir {outputdir} " if outputdir else "")
+        + (f"--json {json} " if json else "")
+        + (f"--petmask {petmask} " if petmask else "")
+        + f"--start {start} "
+        + f"--end {end} "
+        + (f"--fwhm {fwhm} " if fwhm else "")
+        + f"--weight_by {weight_by} "
+        + f"--integration_type {integration_type} "
+        + pet
+    )
+    doc = km.__class__.__doc__
+    derivative_json_dict = {
+        "Description": re.sub(r"\s+", " ", doc) if doc else "",
+        # "Sources": [pet],
+        "ModelName": model_abbr,
+        "ReferenceRegion": refroi if refroi else refmask,
+        "AdditionalModelDetails": (
+            f"Frame weighting by: {weight_by}. "
+            + f"Integration type: {integration_type}.",
+        ),
+        "InputValues": inputvalues,
+        "InputValuesLabels": inputvalueslabels,
+        "InputValuesUnits": inputvaluesunits,
+        "SoftwareName": "dynamicpet",
+        "SoftwareVersion": __version__,
+        "CommandLine": cmd,
+    }
+
+    with open(output_json, "w") as f:
+        json_dump(derivative_json_dict, f, indent=4)
 
 
 def parse_kineticmodel_inputs(
diff --git a/src/dynamicpet/kineticmodel/kineticmodel.py b/src/dynamicpet/kineticmodel/kineticmodel.py
index 0fe6cb3..1650baf 100644
--- a/src/dynamicpet/kineticmodel/kineticmodel.py
+++ b/src/dynamicpet/kineticmodel/kineticmodel.py
@@ -35,6 +35,7 @@ class KineticModel(ABC):
     @abstractmethod
     def get_param_names(cls) -> list[str]:
         """Get names of kinetic model parameters."""
+        # parameter names should contain alphanumeric characters only
         raise NotImplementedError
 
     def __init__(
@@ -110,11 +111,11 @@ def get_parameter(self, param_name: str) -> SpatialImage | NumpyNumberArray:
             else:
                 param_vector: NumpyNumberArray = self.parameters[param_name]
                 return param_vector
-        elif param_name == "bp" and "dvr" in self.parameters:
-            self.parameters[param_name] = self.parameters["dvr"] - 1
+        elif param_name == "BPND" and "DVR" in self.parameters:
+            self.parameters[param_name] = self.parameters["DVR"] - 1
             return self.get_parameter(param_name)
-        elif param_name == "dvr" and "bp" in self.parameters:
-            self.parameters[param_name] = self.parameters["bp"] + 1
+        elif param_name == "DVR" and "BPND" in self.parameters:
+            self.parameters[param_name] = self.parameters["BPND"] + 1
             return self.get_parameter(param_name)
         else:
             raise AttributeError(
diff --git a/src/dynamicpet/kineticmodel/logan.py b/src/dynamicpet/kineticmodel/logan.py
index a631c92..78877a2 100644
--- a/src/dynamicpet/kineticmodel/logan.py
+++ b/src/dynamicpet/kineticmodel/logan.py
@@ -30,7 +30,7 @@ class LRTM(KineticModel):
     @classmethod
     def get_param_names(cls) -> list[str]:
         """Get names of kinetic model parameters."""
-        return ["dvr"]
+        return ["DVR"]
 
     def fit(  # noqa: max-complexity: 12
         self,
@@ -38,7 +38,7 @@ def fit(  # noqa: max-complexity: 12
         integration_type: INTEGRATION_TYPE_OPTS = "trapz",
         weight_by: WEIGHT_OPTS | NumpyNumberArray | None = "frame_duration",
         tstar: float = 0,
-        k2: float | None = None,
+        k2prime: float | None = None,
     ) -> None:
         """Estimate model parameters.
 
@@ -56,7 +56,8 @@ def fit(  # noqa: max-complexity: 12
                   to fit the kinetic model. Elements outside the mask will
                   be set to to 0 in parametric estimate outputs.
             tstar: time beyond which to assume linearity
-            k2: (avg.) effective tissue-to-plasma efflux constant, in unit of 1/min
+            k2prime: (avg.) effective tissue-to-plasma efflux constant in the
+                     reference region, in unit of 1/min
         """
         # get reference TAC as a 1-D vector
         reftac: NumpyNumberArray = self.reftac.dataobj.flatten()[:, np.newaxis]
@@ -82,7 +83,7 @@ def fit(  # noqa: max-complexity: 12
 
         dvr = np.zeros((num_elements, 1))
 
-        if not k2:
+        if not k2prime:
             # TODO
             # Check Eq. 7 assumption (i.e., that tac / reftac is reasonably
             # constant) by calculating R2 etc. for each tac.
@@ -103,13 +104,14 @@ def fit(  # noqa: max-complexity: 12
 
             # ----- Get DVR -----
             # Set up the weighted linear regression model based on Logan et al.:
-            # - use Eq. 6 if k2 is provided
-            # - use Eq. 7 if k2 is not provided
+            # - use Eq. 6 if k2prime is provided
+            # - use Eq. 7 if k2prime is not provided
 
             x = np.column_stack(
                 (
                     np.ones_like(tac_tstar),
-                    (int_reftac_tstar + (reftac_tstar / k2 if k2 else 0)) / tac_tstar,
+                    (int_reftac_tstar + (reftac_tstar / k2prime if k2prime else 0))
+                    / tac_tstar,
                 )
             )
             y = int_tac_tstar / tac_tstar
@@ -123,8 +125,8 @@ def fit(  # noqa: max-complexity: 12
             # distribution volume ratio
             dvr[k] = b[1]
 
-        self.set_parameter("dvr", dvr, mask)
-        # should tstar (and k2?) also be stored?
+        self.set_parameter("DVR", dvr, mask)
+        # should tstar (and k2prime?) also be stored?
 
     def fitted_tacs(self) -> TemporalMatrix | TemporalImage:
         """Get fitted TACs based on estimated model parameters."""
diff --git a/src/dynamicpet/kineticmodel/srtm.py b/src/dynamicpet/kineticmodel/srtm.py
index 269cef4..48ef966 100644
--- a/src/dynamicpet/kineticmodel/srtm.py
+++ b/src/dynamicpet/kineticmodel/srtm.py
@@ -30,7 +30,7 @@ class SRTMLammertsma1996(KineticModel):
     @classmethod
     def get_param_names(cls) -> list[str]:
         """Get names of kinetic model parameters."""
-        return ["bp", "r1", "k2"]
+        return ["BPND", "R1", "k2"]
 
     def fit(
         self,
@@ -56,7 +56,7 @@ def fit(
 
         weights = tacs.get_weights(weight_by)
 
-        bp = np.zeros((num_elements, 1))
+        bp_nd = np.zeros((num_elements, 1))
         r1 = np.zeros((num_elements, 1))
         k2 = np.zeros((num_elements, 1))
         for k in trange(num_elements):
@@ -69,10 +69,10 @@ def fit(
                 sigma=weights,
                 bounds=([0, 0, 0], [15, 10, 1]),
             )
-            bp[k], r1[k], k2[k] = popt
+            bp_nd[k], r1[k], k2[k] = popt
 
-        self.set_parameter("bp", bp, mask)
-        self.set_parameter("r1", r1, mask)
+        self.set_parameter("BPND", bp_nd, mask)
+        self.set_parameter("R1", r1, mask)
         self.set_parameter("k2", k2, mask)
 
     def fitted_tacs(self) -> TemporalMatrix | TemporalImage:
@@ -82,11 +82,11 @@ def fitted_tacs(self) -> TemporalMatrix | TemporalImage:
 
         for i in trange(num_elements):
             idx = np.unravel_index(i, self.tacs.shape[:-1])
-            bp = self.parameters["bp"][*idx]
-            r1 = self.parameters["r1"][*idx]
+            bp_nd = self.parameters["BPND"][*idx]
+            r1 = self.parameters["R1"][*idx]
             k2 = self.parameters["k2"][*idx]
-            if bp or r1 or k2:
-                fitted_tacs_dataobj[*idx, :] = srtm_model(self.reftac, bp, r1, k2)
+            if bp_nd or r1 or k2:
+                fitted_tacs_dataobj[*idx, :] = srtm_model(self.reftac, bp_nd, r1, k2)
 
         if isinstance(self.tacs, TemporalImage):
             img = image_maker(fitted_tacs_dataobj, self.tacs.img)
@@ -100,13 +100,13 @@ def fitted_tacs(self) -> TemporalMatrix | TemporalImage:
 
 
 def srtm_model(
-    reftac: TemporalMatrix, bp: float, r1: float, k2: float
+    reftac: TemporalMatrix, bp_nd: float, r1: float, k2: float
 ) -> NumpyNumberArray:
     """SRTM model to generate a target TAC.
 
     Args:
         reftac: reference TAC
-        bp: binding potential
+        bp_nd: binding potential
         r1: relative radiotracer delivery parameter, R1
         k2: k2
 
@@ -130,7 +130,7 @@ def srtm_model(
         t_upsampled, t, reftac.dataobj.astype("float").flatten()
     )
 
-    k2a = k2 / (1 + bp)
+    k2a = k2 / (1 + bp_nd)
     conv_res_upsampled = (
         convolve(reftac_upsampled, np.exp(-k2a * t_upsampled), mode="full")[
             : len(t_upsampled)
@@ -157,15 +157,14 @@ class SRTMZhou2003(KineticModel):
     def get_param_names(cls) -> list[str]:
         """Get names of kinetic model parameters."""
         return [
-            "dvr",
-            "r1",
+            "DVR",
+            "R1",
             "k2",
             "k2a",
-            "r1_lrsc",
-            "k2_lrsc",
-            "k2a_lrsc",
-            "noise_var_eq_dvr",
-            "noise_var_eq_r1",
+            "R1LRSC",
+            "k2LRSC",
+            "k2aLRSC",
+            "noiseVarEqR1",
         ]
 
     def fit(  # noqa: max-complexity: 12
@@ -277,12 +276,11 @@ def fit(  # noqa: max-complexity: 12
 
             r1[k], k2[k], k2a[k] = b
 
-        self.set_parameter("dvr", dvr, mask)
-        self.set_parameter("r1", r1, mask)
+        self.set_parameter("DVR", dvr, mask)
+        self.set_parameter("R1", r1, mask)
         self.set_parameter("k2", k2, mask)
         self.set_parameter("k2a", k2a, mask)
-        self.set_parameter("noise_var_eq_dvr", noise_var_eq_dvr, mask)
-        self.set_parameter("noise_var_eq_r1", noise_var_eq_r1, mask)
+        self.set_parameter("noiseVarEqR1", noise_var_eq_r1, mask)
 
         if do_smooth:
             smooth_r1_mat, smooth_k2_mat, smooth_k2a_mat, h = self.prep_refine_r1(
@@ -315,9 +313,9 @@ def fit(  # noqa: max-complexity: 12
 
                 r1_lrsc[k], k2_lrsc[k], k2a_lrsc[k] = b
 
-            self.set_parameter("r1_lrsc", r1_lrsc, mask)
-            self.set_parameter("k2_lrsc", k2_lrsc, mask)
-            self.set_parameter("k2a_lrsc", k2a_lrsc, mask)
+            self.set_parameter("R1LRSC", r1_lrsc, mask)
+            self.set_parameter("k2LRSC", k2_lrsc, mask)
+            self.set_parameter("k2aLRSC", k2a_lrsc, mask)
 
     def prep_refine_r1(
         self,
@@ -342,7 +340,7 @@ def prep_refine_r1(
         m = 3
 
         smooth_r1_img: SpatialImage = smooth_image(
-            self.get_parameter("r1"), fwhm  # type: ignore
+            self.get_parameter("R1"), fwhm  # type: ignore
         )
         smooth_k2_img: SpatialImage = smooth_image(
             self.get_parameter("k2"), fwhm  # type: ignore
@@ -351,8 +349,8 @@ def prep_refine_r1(
             self.get_parameter("k2a"), fwhm  # type: ignore
         )
 
-        noise_var_eq_r1_img: SpatialImage = self.get_parameter("noise_var_eq_r1")  # type: ignore
-        r1_img: SpatialImage = self.get_parameter("r1")  # type: ignore
+        noise_var_eq_r1_img: SpatialImage = self.get_parameter("noiseVarEqR1")  # type: ignore
+        r1_img: SpatialImage = self.get_parameter("R1")  # type: ignore
         k2_img: SpatialImage = self.get_parameter("k2")  # type: ignore
         k2a_img: SpatialImage = self.get_parameter("k2a")  # type: ignore
         noise_var_eq_r1_data = noise_var_eq_r1_img.get_fdata()
@@ -409,19 +407,19 @@ def fitted_tacs(self) -> TemporalMatrix | TemporalImage:
         num_elements = self.tacs.num_elements
         fitted_tacs_dataobj = np.empty_like(self.tacs.dataobj)
 
-        use_lrsc = "r1_lrsc" in self.parameters
+        use_lrsc = "R1LRSC" in self.parameters
 
         for i in trange(num_elements):
             idx = np.unravel_index(i, self.tacs.shape[:-1])
-            dvr = self.parameters["dvr"][*idx]
-            r1 = self.parameters["r1"][*idx]
+            dvr = self.parameters["DVR"][*idx]
+            r1 = self.parameters["R1"][*idx]
             k2 = self.parameters["k2"][*idx]
             if dvr or r1 or k2:
-                bp = dvr - 1
+                bp_nd = dvr - 1
                 if use_lrsc:
-                    r1 = self.parameters["r1_lrsc"][*idx]
-                    k2 = self.parameters["k2_lrsc"][*idx]
-                fitted_tacs_dataobj[*idx, :] = srtm_model(self.reftac, bp, r1, k2)
+                    r1 = self.parameters["R1LRSC"][*idx]
+                    k2 = self.parameters["k2LRSC"][*idx]
+                fitted_tacs_dataobj[*idx, :] = srtm_model(self.reftac, bp_nd, r1, k2)
 
         if isinstance(self.tacs, TemporalImage):
             img = image_maker(fitted_tacs_dataobj, self.tacs.img)
diff --git a/src/dynamicpet/kineticmodel/suvr.py b/src/dynamicpet/kineticmodel/suvr.py
index d3452c7..ea9fbcf 100644
--- a/src/dynamicpet/kineticmodel/suvr.py
+++ b/src/dynamicpet/kineticmodel/suvr.py
@@ -20,7 +20,7 @@ class SUVR(KineticModel):
     @classmethod
     def get_param_names(cls) -> list[str]:
         """Get names of kinetic model parameters."""
-        return ["suvr"]
+        return ["SUVR"]
 
     def fit(self, mask: NumpyNumberArray | None = None) -> None:
         """Calculate SUVR.
@@ -45,7 +45,7 @@ def fit(self, mask: NumpyNumberArray | None = None) -> None:
                                       frame_duration=frame_duration)
             >>> km = SUVR(reftac, tacs)
             >>> km.fit()
-            >>> km.get_parameter('suvr')
+            >>> km.get_parameter('SUVR')
             array([1.5, 3. ])
         """
         tacs: TemporalMatrix = self.tacs.timeseries_in_mask(mask)
@@ -58,7 +58,7 @@ def fit(self, mask: NumpyNumberArray | None = None) -> None:
         )
         suvr = numerator / denominator
 
-        self.set_parameter("suvr", suvr, mask)
+        self.set_parameter("SUVR", suvr, mask)
 
     def fitted_tacs(self) -> TemporalMatrix | TemporalImage:
         """Get fitted TACs based on estimated model parameters."""
diff --git a/src/dynamicpet/petbids/petbidsimage.py b/src/dynamicpet/petbids/petbidsimage.py
index ee5d40c..9f36c4a 100644
--- a/src/dynamicpet/petbids/petbidsimage.py
+++ b/src/dynamicpet/petbids/petbidsimage.py
@@ -126,25 +126,29 @@ def decay_uncorrect(self) -> "PETBIDSImage":
     def to_filename(
         self,
         filename: str | PathLike[str],
+        save_json: bool = False,
         anchor: Literal["InjectionStart", "ScanStart"] = "InjectionStart",
     ) -> None:
         """Save to file.
 
         Args:
             filename: file name for the PET image output
+            save_json: whether the PET-BIDS json side car should be saved
             anchor: time anchor. The corresponding tag in the PET-BIDS json will
                     be set to zero (with appropriate offsets applied to other
                     tags).
         """
-        self.set_timezero(anchor)
         self.img.to_filename(filename)
 
-        fbase, fext = op.splitext(filename)
-        if fext == ".gz":
-            fbase = op.splitext(fbase)[0]
-        jsonfilename = fbase + ".json"
+        if save_json:
+            self.set_timezero(anchor)
 
-        write_json(self.json_dict, jsonfilename)
+            fbase, fext = op.splitext(filename)
+            if fext == ".gz":
+                fbase = op.splitext(fbase)[0]
+            jsonfilename = fbase + ".json"
+
+            write_json(self.json_dict, jsonfilename)
 
 
 def load(
diff --git a/src/dynamicpet/petbids/petbidsjson.py b/src/dynamicpet/petbids/petbidsjson.py
index 12cf0c4..5ac072c 100644
--- a/src/dynamicpet/petbids/petbidsjson.py
+++ b/src/dynamicpet/petbids/petbidsjson.py
@@ -289,4 +289,4 @@ def write_json(json_dict: PetBidsJson, filename: str | PathLike[str]) -> None:
         filename: file name for the output
     """
     with open(filename, "w") as f:
-        json_dump(json_dict, f)
+        json_dump(json_dict, f, indent=4)
diff --git a/src/dynamicpet/petbids/petbidsmatrix.py b/src/dynamicpet/petbids/petbidsmatrix.py
index 28868d8..6686dcf 100644
--- a/src/dynamicpet/petbids/petbidsmatrix.py
+++ b/src/dynamicpet/petbids/petbidsmatrix.py
@@ -131,12 +131,14 @@ def decay_uncorrect(self) -> "PETBIDSMatrix":
     def to_filename(
         self,
         filename: str | PathLike[str],
+        save_json: bool = False,
         anchor: Literal["InjectionStart", "ScanStart"] = "InjectionStart",
     ) -> None:
         """Save to file.
 
         Args:
             filename: file name for the tabular TAC tsv output
+            save_json: whether the PET-BIDS json side car should be saved
             anchor: time anchor. The corresponding tag in the PET-BIDS json will
                     be set to zero (with appropriate offsets applied to other
                     tags).
@@ -144,19 +146,20 @@ def to_filename(
         Raises:
             ValueError: file is not a tsv file
         """
-        self.set_timezero(anchor)
-        fbase, fext = op.splitext(filename)
-        if fext != ".tsv":
-            raise ValueError("output file must be a tsv file")
-        jsonfilename = fbase + ".json"
-
         with open(filename, "w") as f:
             tsvwriter = csv.writer(f, delimiter="\t")
             tsvwriter.writerow(self.elem_names)
             for row in self.dataobj.T:
                 tsvwriter.writerow(row)
 
-        write_json(self.json_dict, jsonfilename)
+        if save_json:
+            self.set_timezero(anchor)
+
+            fbase, fext = op.splitext(filename)
+            if fext != ".tsv":
+                raise ValueError("output file must be a tsv file")
+            jsonfilename = fbase + ".json"
+            write_json(self.json_dict, jsonfilename)
 
 
 def load(
diff --git a/tests/test_kineticmodel.py b/tests/test_kineticmodel.py
index 1221f8d..f0a3d61 100644
--- a/tests/test_kineticmodel.py
+++ b/tests/test_kineticmodel.py
@@ -75,7 +75,7 @@ def test_suvr_tm(reftac: TemporalMatrix, tacs: TemporalMatrix) -> None:
     """Test SUVR calculation using TemporalMatrix."""
     km = SUVR(reftac, tacs)
     km.fit()
-    suvr: NumpyRealNumberArray = km.get_parameter("suvr")  # type: ignore
+    suvr: NumpyRealNumberArray = km.get_parameter("SUVR")  # type: ignore
     assert suvr.shape == (2,)
     assert np.all(suvr == [2, 3])
 
@@ -87,7 +87,7 @@ def test_suvr_ti(reftac: TemporalMatrix, tacs_img: TemporalImage) -> None:
 
     fitted_tacs = km.fitted_tacs()
 
-    suvr_img: SpatialImage = km.get_parameter("suvr")  # type: ignore
+    suvr_img: SpatialImage = km.get_parameter("SUVR")  # type: ignore
     assert suvr_img.shape == (1, 1, 2)
     assert np.all(suvr_img.get_fdata() == np.array([[[2, 3]]]))
     assert np.all(fitted_tacs.dataobj == tacs_img.dataobj)
@@ -97,13 +97,13 @@ def test_srtm_zhou2003_ti(reftac: TemporalMatrix, tacs_img: TemporalImage) -> No
     """Test SRTM Zhou 2003 using TemporalImage."""
     km = SRTMZhou2003(reftac, tacs_img)
     km.fit(integration_type="trapz")
-    dvr_img: SpatialImage = km.get_parameter("dvr")  # type: ignore
+    dvr_img: SpatialImage = km.get_parameter("DVR")  # type: ignore
 
     assert dvr_img.shape == (1, 1, 2)
 
-    bp_img: SpatialImage = km.get_parameter("bp")  # type: ignore
+    bp_nd_img: SpatialImage = km.get_parameter("BPND")  # type: ignore
 
-    assert np.allclose(dvr_img.get_fdata(), bp_img.get_fdata() + 1)
+    assert np.allclose(dvr_img.get_fdata(), bp_nd_img.get_fdata() + 1)
 
 
 def test_logan_tm(reftac: TemporalMatrix, tacs_img: TemporalImage) -> None:
@@ -111,5 +111,5 @@ def test_logan_tm(reftac: TemporalMatrix, tacs_img: TemporalImage) -> None:
     km = LRTM(reftac, tacs_img)
     km.fit(integration_type="trapz")
 
-    dvr_img: SpatialImage = km.get_parameter("dvr")  # type: ignore
+    dvr_img: SpatialImage = km.get_parameter("DVR")  # type: ignore
     assert dvr_img.shape == (1, 1, 2)
diff --git a/tests/test_kineticmodel_simulated.py b/tests/test_kineticmodel_simulated.py
index 75cb1ef..deb3321 100644
--- a/tests/test_kineticmodel_simulated.py
+++ b/tests/test_kineticmodel_simulated.py
@@ -108,36 +108,36 @@ def frame_duration() -> NDArray[np.double]:
 def srtm_gen_ode_model(
     y: NumpyRealNumberArray,
     t: NumpyRealNumberArray,
-    bp_val: NumpyRealNumberArray,
+    bp_nd_val: NumpyRealNumberArray,
     r1_val: NumpyRealNumberArray,
     plasma_rate: float = -0.03,
-    kref1: float = 1.0,
-    kref2: float = 0.3,
+    k1prime: float = 1.0,
+    k2prime: float = 0.3,
 ) -> list[Any]:
     """Generative ODE model describing SRTM.
 
     Args:
         y: concentrations
         t: time
-        bp_val: binding potential
+        bp_nd_val: binding potential
         r1_val: relative radiotracer delivery
         plasma_rate: plasma rate
-        kref1: k1 for reference region
-        kref2: k2 for reference region
+        k1prime: k1 for reference region
+        k2prime: k2 for reference region
 
     Returns:
         time derivatives of concentrations
     """
-    k1 = r1_val * kref1
+    k1 = r1_val * k1prime
 
-    dvr = 1 + bp_val
-    k2a = kref2 * r1_val / dvr
+    dvr = 1 + bp_nd_val
+    k2a = k2prime * r1_val / dvr
 
     cp, ct, cr = y
 
     dcpdt = plasma_rate * cp
     dctdt = k1 * cp - k2a * ct
-    dcrdt = kref1 * cp - kref2 * cr
+    dcrdt = k1prime * cp - k2prime * cr
 
     dydt = [dcpdt, dctdt, dcrdt]
 
@@ -147,7 +147,7 @@ def srtm_gen_ode_model(
 def get_tacs_and_reftac_dataobj(
     frame_start: NDArray[np.double],
     frame_duration: NDArray[np.double],
-    bp: float,
+    bp_nd: float,
     r1: float,
 ) -> tuple[NDArray[np.double], NDArray[np.double]]:
     """Get digitized TACs."""
@@ -160,7 +160,7 @@ def get_tacs_and_reftac_dataobj(
     num_ode_pts = 500
 
     t_ode = np.linspace(frame_start[0], frame_end[-1], num_ode_pts)
-    y = odeint(srtm_gen_ode_model, y0, t_ode, args=(bp, r1))
+    y = odeint(srtm_gen_ode_model, y0, t_ode, args=(bp_nd, r1))
 
     # Digitize this curve
     cref = np.zeros(len(frame_start))
@@ -178,10 +178,10 @@ def test_srtmzhou2003_tm(
     frame_duration: NDArray[np.double],
 ) -> None:
     """Test SRTM Zhou 2003 calculation using TemporalMatrix."""
-    bp_true = 1.5
+    bp_nd_true = 1.5
     r1_true = 1.2
     ct, cref = get_tacs_and_reftac_dataobj(
-        frame_start, frame_duration, bp_true, r1_true
+        frame_start, frame_duration, bp_nd_true, r1_true
     )
 
     tac = TemporalMatrix(ct, frame_start, frame_duration)
@@ -194,11 +194,11 @@ def test_srtmzhou2003_tm(
 
     fitted_tacs = km.fitted_tacs()
 
-    bp: NumpyRealNumberArray = km.get_parameter("bp")  # type: ignore
-    r1: NumpyRealNumberArray = km.get_parameter("r1")  # type: ignore
+    bp_nd: NumpyRealNumberArray = km.get_parameter("BPND")  # type: ignore
+    r1: NumpyRealNumberArray = km.get_parameter("R1")  # type: ignore
 
     relative_tol = 0.007  # .007 means that 0.7% error is tolerated
-    assert np.allclose(bp, bp_true, rtol=relative_tol)
+    assert np.allclose(bp_nd, bp_nd_true, rtol=relative_tol)
     assert np.allclose(r1, r1_true, rtol=relative_tol)
     assert np.allclose(fitted_tacs.dataobj, tac.dataobj, rtol=0.01)
 
@@ -208,10 +208,10 @@ def test_srtmzhou2003_ti(
     frame_duration: NDArray[np.double],
 ) -> None:
     """Test SRTM Zhou 2003 calculation using TemporalImage."""
-    bp_true = 1.5
+    bp_nd_true = 1.5
     r1_true = 1.2
     ct, cref = get_tacs_and_reftac_dataobj(
-        frame_start, frame_duration, bp_true, r1_true
+        frame_start, frame_duration, bp_nd_true, r1_true
     )
 
     dims = (10, 11, 12, len(frame_start))
@@ -232,12 +232,12 @@ def test_srtmzhou2003_ti(
 
     fitted_tacs = km.fitted_tacs()
 
-    bp: SpatialImage = km.get_parameter("bp")  # type: ignore
-    r1: SpatialImage = km.get_parameter("r1")  # type: ignore
-    r1_lrsc: SpatialImage = km.get_parameter("r1_lrsc")  # type: ignore
+    bp_nd: SpatialImage = km.get_parameter("BPND")  # type: ignore
+    r1: SpatialImage = km.get_parameter("R1")  # type: ignore
+    r1_lrsc: SpatialImage = km.get_parameter("R1LRSC")  # type: ignore
 
     relative_tol = 0.007  # .007 means that 0.7% error is tolerated
-    assert np.allclose(bp.get_fdata(), bp_true, rtol=relative_tol)
+    assert np.allclose(bp_nd.get_fdata(), bp_nd_true, rtol=relative_tol)
     assert np.allclose(r1.get_fdata(), r1_true, rtol=relative_tol)
     assert np.allclose(r1_lrsc.get_fdata(), r1_true, rtol=relative_tol)
     assert np.allclose(fitted_tacs.dataobj, tac.dataobj, rtol=0.01)
@@ -248,10 +248,10 @@ def test_srtmlammertsma1996_tm(
     frame_duration: NDArray[np.double],
 ) -> None:
     """Test SRTM Lammertsma 1996 calculation using TemporalMatrix."""
-    bp_true = 1.5
+    bp_nd_true = 1.5
     r1_true = 1.2
     ct, cref = get_tacs_and_reftac_dataobj(
-        frame_start, frame_duration, bp_true, r1_true
+        frame_start, frame_duration, bp_nd_true, r1_true
     )
 
     tac = TemporalMatrix(ct, frame_start, frame_duration)
@@ -262,11 +262,11 @@ def test_srtmlammertsma1996_tm(
 
     fitted_tacs = km.fitted_tacs()
 
-    bp: NumpyRealNumberArray = km.get_parameter("bp")  # type: ignore
-    r1: NumpyRealNumberArray = km.get_parameter("r1")  # type: ignore
+    bp_nd: NumpyRealNumberArray = km.get_parameter("BPND")  # type: ignore
+    r1: NumpyRealNumberArray = km.get_parameter("R1")  # type: ignore
 
     relative_tol = 0.02  # .02 means that 2% error is tolerated
-    assert np.allclose(bp, bp_true, rtol=relative_tol)
+    assert np.allclose(bp_nd, bp_nd_true, rtol=relative_tol)
     assert np.allclose(r1, r1_true, rtol=relative_tol)
     assert np.allclose(fitted_tacs.dataobj, tac.dataobj, rtol=relative_tol)
 
@@ -276,10 +276,10 @@ def test_lrtm_tm(
     frame_duration: NDArray[np.double],
 ) -> None:
     """Test Logan 1996 calculation using TemporalMatrix."""
-    bp_true = 1.5
+    bp_nd_true = 1.5
     r1_true = 1.2
     ct, cref = get_tacs_and_reftac_dataobj(
-        frame_start, frame_duration, bp_true, r1_true
+        frame_start, frame_duration, bp_nd_true, r1_true
     )
 
     tac = TemporalMatrix(ct, frame_start, frame_duration)
@@ -290,7 +290,7 @@ def test_lrtm_tm(
     # to fit to obtain the best possible recovery of true values
     km.fit(integration_type="trapz")
 
-    bp: NumpyRealNumberArray = km.get_parameter("bp")  # type: ignore
+    bp_nd: NumpyRealNumberArray = km.get_parameter("BPND")  # type: ignore
 
     relative_tol = 0.02  # .02 means that 2% error is tolerated
-    assert np.allclose(bp, bp_true, rtol=relative_tol)
+    assert np.allclose(bp_nd, bp_nd_true, rtol=relative_tol)
diff --git a/tests/test_main.py b/tests/test_main.py
index df1de48..0a90a93 100644
--- a/tests/test_main.py
+++ b/tests/test_main.py
@@ -1,7 +1,6 @@
 """Test cases for __main__.py."""
 
 import csv
-import os
 from pathlib import Path
 
 import numpy as np
@@ -32,9 +31,9 @@ def images() -> dict[str, Path]:
     outdir = Path(__file__).parent / "test_data"
     outdir.mkdir(exist_ok=True)
 
-    rois_fname = outdir / "displacementROI_dseg.nii.gz"
-    petjson_fname = outdir / "pet.json"
-    pet_fname = outdir / "pet.nii"
+    rois_fname = outdir / "sub-000101_ses-baseline_label-displacementROI_dseg.nii.gz"
+    petjson_fname = outdir / "sub-000101_ses-baseline_pet.json"
+    pet_fname = outdir / "sub-000101_ses-baseline_pet.nii"
 
     baseurl = "https://s3.amazonaws.com/openneuro.org/ds001705/sub-000101/ses-baseline/"
     roisurl = (
@@ -72,9 +71,9 @@ def images() -> dict[str, Path]:
                 for chunk in r.iter_content(chunk_size=8192):
                     f.write(chunk)
 
-    petmask_fname = outdir / "petmask.nii.gz"
-    refmask_fname = outdir / "refmask.nii.gz"
-    rois_csv_fname = outdir / "rois.tsv"
+    petmask_fname = outdir / "sub-000101_ses-baseline_desc-pet_mask.nii.gz"
+    refmask_fname = outdir / "sub-000101_ses-baseline_desc-ref_mask.nii.gz"
+    rois_csv_fname = outdir / "sub-000101_ses-baseline_tacs.tsv"
 
     if petmask_fname.exists() and refmask_fname.exists() and rois_csv_fname.exists():
         pass
@@ -142,8 +141,8 @@ def test_denoise_hyprlr(images: dict[str, Path]) -> None:
     )
 
     assert result.exit_code == 0
-    assert os.path.isfile(outputdir / "pet_hyprlr.nii")
-    assert os.path.isfile(outputdir / "pet_hyprlr.json")
+    assert (outputdir / "sub-000101_ses-baseline_desc-hyprlr_pet.nii").is_file()
+    assert (outputdir / "sub-000101_ses-baseline_desc-hyprlr_pet.json").is_file()
 
 
 def test_denoise_nesma(images: dict[str, Path]) -> None:
@@ -176,8 +175,8 @@ def test_denoise_nesma(images: dict[str, Path]) -> None:
     )
 
     assert result.exit_code == 0
-    assert os.path.isfile(outputdir / "pet_nesma.nii")
-    assert os.path.isfile(outputdir / "pet_nesma.json")
+    assert (outputdir / "sub-000101_ses-baseline_desc-nesma_pet.nii").is_file()
+    assert (outputdir / "sub-000101_ses-baseline_desc-nesma_pet.json").is_file()
 
 
 def test_kineticmodel_suvr(images: dict[str, Path]) -> None:
@@ -210,7 +209,7 @@ def test_kineticmodel_suvr(images: dict[str, Path]) -> None:
         catch_exceptions=False,
     )
     assert result.exit_code == 0
-    assert (outputdir / "rois_km-suvr.tsv").is_file()
+    assert (outputdir / "sub-000101_ses-baseline_model-SUVR_kinpar.tsv").is_file()
 
     # next, test with nifti images
     pet_fname = images["pet_fname"]
@@ -239,21 +238,27 @@ def test_kineticmodel_suvr(images: dict[str, Path]) -> None:
         catch_exceptions=False,
     )
     assert result.exit_code == 0
-    assert os.path.isfile(outputdir / "pet_km-suvr_kp-suvr.nii")
+    assert (
+        outputdir / "sub-000101_ses-baseline_model-SUVR_meas-SUVR_mimap.nii"
+    ).is_file()
+    assert (outputdir / "sub-000101_ses-baseline_model-SUVR_mimap.json").is_file()
 
     # finally, make sure that the two methods give the same results
 
     # skip the header (column names) and the first column (ROI names)
     rois_km_suvr = np.genfromtxt(
-        outputdir / "rois_km-suvr.tsv", delimiter="\t", skip_header=1, usecols=(1,)
+        outputdir / "sub-000101_ses-baseline_model-SUVR_kinpar.tsv",
+        delimiter="\t",
+        skip_header=1,
+        usecols=(1,),
     )
 
     # read in rois_fname
     rois_img: SpatialImage = nib_load(images["rois_fname"])  # type: ignore
 
-    # calculate the average SUVR (pet_km-suvr_kp-suvr.nii) in each ROI in rois_img
+    # calculate avg SUVR (pet_model-SUVR_meas-SUVR_mimap.nii) per ROI in rois_img
     pet_km_suvr_img: SpatialImage = nib_load(  # type: ignore
-        outputdir / "pet_km-suvr_kp-suvr.nii"
+        outputdir / "sub-000101_ses-baseline_model-SUVR_meas-SUVR_mimap.nii"
     )
     pet_km_suvr = pet_km_suvr_img.get_fdata()
     num_rois = len(np.unique(rois_img.get_fdata())) - 1
@@ -298,7 +303,8 @@ def test_kineticmodel_srtmzhou2003(images: dict[str, Path]) -> None:
         catch_exceptions=False,
     )
     assert result.exit_code == 0
-    assert os.path.isfile(outputdir / ("rois_km-" + model + ".tsv"))
+    assert (outputdir / "sub-000101_ses-baseline_model-SRTM_kinpar.tsv").is_file()
+    assert (outputdir / "sub-000101_ses-baseline_model-SRTM_kinpar.json").is_file()
 
     # next, test with nifti images
     pet_fname = images["pet_fname"]
@@ -327,5 +333,10 @@ def test_kineticmodel_srtmzhou2003(images: dict[str, Path]) -> None:
         catch_exceptions=False,
     )
     assert result.exit_code == 0
-    assert os.path.isfile(outputdir / ("pet_km-" + model + "_kp-dvr.nii"))
-    assert os.path.isfile(outputdir / ("pet_km-" + model + "_kp-r1.nii"))
+    assert (
+        outputdir / "sub-000101_ses-baseline_model-SRTM_meas-DVR_mimap.nii"
+    ).is_file()
+    assert (
+        outputdir / "sub-000101_ses-baseline_model-SRTM_meas-R1_mimap.nii"
+    ).is_file()
+    assert (outputdir / "sub-000101_ses-baseline_model-SRTM_mimap.json").is_file()
diff --git a/tests/test_petbidsimage.py b/tests/test_petbidsimage.py
index 72b2fec..8d52072 100644
--- a/tests/test_petbidsimage.py
+++ b/tests/test_petbidsimage.py
@@ -324,7 +324,7 @@ def test_hypr_lr(ti: PETBIDSImage) -> None:
 def test_file_io(ti: PETBIDSImage, tmp_path: Path) -> None:
     """Test writing to file and reading it back."""
     fname = tmp_path / "test.nii.gz"
-    ti.to_filename(fname)
+    ti.to_filename(fname, save_json=True)
     ti2 = load(fname)
 
     assert np.allclose(ti.frame_start, ti2.frame_start)
diff --git a/tests/test_petbidsmatrix.py b/tests/test_petbidsmatrix.py
index c63af10..a6a3954 100644
--- a/tests/test_petbidsmatrix.py
+++ b/tests/test_petbidsmatrix.py
@@ -52,7 +52,7 @@ def test_extract(pm: PETBIDSMatrix) -> None:
 def test_file_io(pm: PETBIDSMatrix, tmp_path: Path) -> None:
     """Test writing to file and reading it back."""
     fname = tmp_path / "test.tsv"
-    pm.to_filename(fname)
+    pm.to_filename(fname, save_json=True)
     pm2 = load(fname)
 
     assert np.allclose(pm.frame_start, pm2.frame_start)