diff --git a/CHANGELOG.md b/CHANGELOG.md
index 96d963e2b8..fde17e0df2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,6 +15,7 @@
 - Renamed "electrode diffusivity" to "particle diffusivity" as a non-breaking change with a deprecation warning ([#3624](https://github.com/pybamm-team/PyBaMM/pull/3624))
 - Add support for BPX version 0.4.0 which allows for blended electrodes and user-defined parameters in BPX([#3414](https://github.com/pybamm-team/PyBaMM/pull/3414))
 - Added `by_submodel` feature in `print_parameter_info` method to allow users to print parameters and types of submodels in a tabular and readable format ([#3628](https://github.com/pybamm-team/PyBaMM/pull/3628))
+- Added `WyciskOpenCircuitPotential` for differential capacity hysteresis state open-circuit potential submodel ([#3593](https://github.com/pybamm-team/PyBaMM/pull/3593))
 
 ## Bug Fixes
 
diff --git a/docs/source/api/models/submodels/interface/open_circuit_potential/index.rst b/docs/source/api/models/submodels/interface/open_circuit_potential/index.rst
index fc664adf2b..3a20bccb17 100644
--- a/docs/source/api/models/submodels/interface/open_circuit_potential/index.rst
+++ b/docs/source/api/models/submodels/interface/open_circuit_potential/index.rst
@@ -7,3 +7,4 @@ Open-circuit potential models
   current_sigmoid_ocp
   single_ocp
   msmr_ocp
+  wycisk_ocp
diff --git a/docs/source/api/models/submodels/interface/open_circuit_potential/wycisk_ocp.rst b/docs/source/api/models/submodels/interface/open_circuit_potential/wycisk_ocp.rst
new file mode 100644
index 0000000000..6952a2ee44
--- /dev/null
+++ b/docs/source/api/models/submodels/interface/open_circuit_potential/wycisk_ocp.rst
@@ -0,0 +1,5 @@
+Wycisk Open Circuit Potential
+=============================
+
+.. autoclass:: pybamm.open_circuit_potential.WyciskOpenCircuitPotential
+    :members:
diff --git a/docs/source/examples/index.rst b/docs/source/examples/index.rst
index 8501e9cba8..a4340adf8d 100644
--- a/docs/source/examples/index.rst
+++ b/docs/source/examples/index.rst
@@ -50,6 +50,7 @@ The notebooks are organised into subfolders, and can be viewed in the galleries
     notebooks/models/compare-particle-diffusion-models.ipynb
     notebooks/models/composite_particle.ipynb
     notebooks/models/coupled-degradation.ipynb
+    notebooks/models/differential-capacity-hysteresis-state.ipynb
     notebooks/models/DFN-with-particle-size-distributions.ipynb
     notebooks/models/DFN.ipynb
     notebooks/models/electrode-state-of-health.ipynb
diff --git a/docs/source/examples/notebooks/models/differential-capacity-hysteresis-state.ipynb b/docs/source/examples/notebooks/models/differential-capacity-hysteresis-state.ipynb
new file mode 100644
index 0000000000..3a3bb81a53
--- /dev/null
+++ b/docs/source/examples/notebooks/models/differential-capacity-hysteresis-state.ipynb
@@ -0,0 +1,269 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Differential Capacity Hysteresis State model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Note: you may need to restart the kernel to use updated packages.\n"
+     ]
+    }
+   ],
+   "source": [
+    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
+    "import pybamm"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Model Equations\n",
+    "\n",
+    "Herein the model equations for the Differential Capacity Hysteresis State open-circuit potential model are outlined, as described in Wycisk (2022).\n",
+    "\n",
+    "### Hysteresis State Variable\n",
+    "\n",
+    "This approach utilizes a state variable to represent the degree of hysteresis at a given time and stoichiometry, $h(z,t)$. The hysteresis is treated separately from the open-circuit potential, where the potential of the electrode is written as\n",
+    "\n",
+    "$$ U = U_{avg}^0(z) + H(z) \\cdot h(z,t) - \\eta $$\n",
+    "\n",
+    "Where $H(z)$ is a function representing the hysteresis as a function of stoichiometry, $z$, and where $\\eta$ represents the sum of the overpotentials. $U_{avg}^0(z)$ is simply the average of the delithiation and lithiation open-circuit potential branches. $H(z)$ can be determined by finding the half-difference value between the lithiation and delithiation branches across the entire stoichiometry range. The state variable $h(z,t)$ is both stoichiometry and time-dependant, and spans between the range of -1 and 1. The hysteresis state variable $h(z,t)$ can be expressed in differential form with respect to time as\n",
+    "\n",
+    "$$ \\frac{dh(z,t)}{dt} = \\left(\\frac{k(z) \\cdot I(t)}{Q_{cell}}\\right)\\left(1-\\text{sgn}\\left(\\frac{dz(t)}{dt}\\right) h(z,t)\\right) $$\n",
+    "\n",
+    "where $ k(z) $ is expressed as \n",
+    "\n",
+    "$$ k(z) = K \\cdot \\frac{1}{\\left(C_{diff}\\left(z\\right)\\right)^{x}} $$\n",
+    "\n",
+    "And where $C_{diff}(z)$ is the differential capacity with respect to potential, expressed as \n",
+    "\n",
+    "$$ C_{diff}(z) = \\frac{dQ}{dU_{avg}^0(z)} $$\n",
+    "\n",
+    "Here, $Q$ is the capacity of the phase or active material experiencing the voltage hysteresis. The remaining parameters are $K$ and $x$ which are both fitting parameters that affect the response of the hysteresis state decay when passing charge in either direction.\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Comparing the DCHS and Current-Sigmoid model approaches\n",
+    "\n",
+    "The behavior of the DCHS model is different than the current-sigmoid model approach for open-circuit potential in systems with hysteresis. Where the current-sigmoid model switches between hysteresis states simply based on the instantaneous current, the DCHS model switches based on the amount of charge passed through the active material phase while also relying on the previous hysteresis state. To assess this differentiated performance, we will compare it to the current-sigmoid model by adapting the Chen2020_composite parameter set.\n",
+    "\n",
+    "First we generate the model, and specify the open-circuit potential methods for the negative and positive electrodes. To maintain consistency with the parameter set, two phases for the negative electrode will be defined."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model_DCHS = pybamm.lithium_ion.DFN(\n",
+    "    {\n",
+    "        \"open-circuit potential\": ((\"single\", \"Wycisk\"), \"single\"),\n",
+    "        \"particle phases\": (\"2\", \"1\"),\n",
+    "    }\n",
+    ")\n",
+    "\n",
+    "model_current_sigmoid = pybamm.lithium_ion.DFN(\n",
+    "    {\n",
+    "        \"open-circuit potential\": ((\"single\", \"current sigmoid\"), \"single\"),\n",
+    "        \"particle phases\": (\"2\", \"1\"),\n",
+    "    }\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, lets define the modifications to the parameter set"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "parameters_DCHS = pybamm.ParameterValues(\"Chen2020_composite\")\n",
+    "parameters_current_sigmoid = pybamm.ParameterValues(\"Chen2020_composite\")\n",
+    "\n",
+    "\n",
+    "# get the lithiation and delithiation functions\n",
+    "lithiation_ocp = parameters_DCHS[\"Secondary: Negative electrode lithiation OCP [V]\"]\n",
+    "delithiation_ocp = parameters_DCHS[\"Secondary: Negative electrode delithiation OCP [V]\"]\n",
+    "\n",
+    "\n",
+    "# define an additional OCP function\n",
+    "def ocp_avg(sto):\n",
+    "    return (lithiation_ocp(sto) + delithiation_ocp(sto)) / 2\n",
+    "\n",
+    "\n",
+    "# add additional parameters\n",
+    "parameters_DCHS.update(\n",
+    "    {\n",
+    "        \"Secondary: Negative electrode OCP [V]\": ocp_avg,\n",
+    "    },\n",
+    "    check_already_exists=False,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, we need to add the additional parameters required by the model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "parameters_DCHS.update(\n",
+    "    {\n",
+    "        \"Secondary: Negative particle hysteresis decay rate\": 0.005,\n",
+    "        \"Secondary: Negative particle hysteresis switching factor\": 10,\n",
+    "    },\n",
+    "    check_already_exists=False,\n",
+    ")\n",
+    "\n",
+    "\n",
+    "experiment = pybamm.Experiment(\n",
+    "    [\n",
+    "        (\"Discharge at 1 C for 1 hour or until 2.5 V\", \"Rest for 15 minutes\"),\n",
+    "        (\n",
+    "            \"Charge at 1C until 4.2 V\",\n",
+    "            \"Hold at 4.2 V until 0.05 C\",\n",
+    "            \"Rest for 15 minutes\",\n",
+    "        ),\n",
+    "    ]\n",
+    ")\n",
+    "\n",
+    "\n",
+    "simulation_dchs = pybamm.Simulation(\n",
+    "    model_DCHS, experiment=experiment, parameter_values=parameters_DCHS\n",
+    ")\n",
+    "solution_dchs = simulation_dchs.solve(calc_esoh=False)\n",
+    "\n",
+    "simulation_current_sigmoid = pybamm.Simulation(\n",
+    "    model_current_sigmoid,\n",
+    "    experiment=experiment,\n",
+    "    parameter_values=parameters_current_sigmoid,\n",
+    ")\n",
+    "\n",
+    "solution_current_sigmoid = simulation_current_sigmoid.solve(calc_esoh=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now plotting the results and the hysteresis state "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "e6677ed985c14dd8941223b20650f6fd",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=3.1492654802910014, step=0.03149265480291001…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "output_variables = [\n",
+    "    \"X-averaged negative electrode secondary hysteresis state\",\n",
+    "    \"Negative electrode secondary open-circuit potential [V]\",\n",
+    "    \"Negative electrode secondary stoichiometry\",\n",
+    "    \"Terminal voltage [V]\",\n",
+    "    \"X-averaged negative electrode secondary open-circuit potential [V]\",\n",
+    "]\n",
+    "\n",
+    "pybamm.QuickPlot(solution_dchs, output_variables=output_variables).dynamic_plot()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "51ea98d2812c4afd97b1b9c33ee95eef",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=3.1492654802910014, step=0.03149265480291001…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "output_variables = [\n",
+    "    \"Terminal voltage [V]\",\n",
+    "    \"Current [A]\",\n",
+    "    \"Negative electrode secondary open-circuit potential [V]\",\n",
+    "]\n",
+    "pybamm.QuickPlot(\n",
+    "    [solution_current_sigmoid, solution_dchs],\n",
+    "    labels=[\"Current sigmoid\", \"DCHS\"],\n",
+    "    output_variables=output_variables,\n",
+    ").dynamic_plot()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "venv",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/pybamm/CITATIONS.bib b/pybamm/CITATIONS.bib
index d2a6643f69..d19093a323 100644
--- a/pybamm/CITATIONS.bib
+++ b/pybamm/CITATIONS.bib
@@ -691,3 +691,16 @@ @article{landesfeind2019temperature
   year={2019},
   publisher={The Electrochemical Society}
 }
+
+@article{Wycisk2022,
+title = {Modified Plett-model for modeling voltage hysteresis in lithium-ion cells},
+journal = {Journal of Energy Storage},
+volume = {52},
+pages = {105016},
+year = {2022},
+issn = {2352-152X},
+doi = {https://doi.org/10.1016/j.est.2022.105016},
+url = {https://www.sciencedirect.com/science/article/pii/S2352152X22010192},
+author = {Dominik Wycisk and Marc Oldenburger and Marc Gerry Stoye and Toni Mrkonjic and Arnulf Latz},
+keywords = {Lithium-ion battery, Voltage hysteresis, Plett-model, Silicon–graphite anode},
+}
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index e12f211abb..ed2a490783 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -102,7 +102,7 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 reactions.
             * "open-circuit potential" : str
                 Sets the model for the open circuit potential. Can be "single"
-                (default), "current sigmoid", or "MSMR". If "MSMR" then the "particle"
+                (default), "current sigmoid", "Wycisk", or "MSMR". If "MSMR" then the "particle"
                 option must also be "MSMR". A 2-tuple can be provided for different
                 behaviour in negative and positive electrodes.
             * "operating mode" : str
@@ -263,7 +263,7 @@ def __init__(self, extra_options):
                 "stress and reaction-driven",
             ],
             "number of MSMR reactions": ["none"],
-            "open-circuit potential": ["single", "current sigmoid", "MSMR"],
+            "open-circuit potential": ["single", "current sigmoid", "MSMR", "Wycisk"],
             "operating mode": [
                 "current",
                 "voltage",
diff --git a/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py b/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
index 879b8effe1..479e8203ed 100644
--- a/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
+++ b/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
@@ -248,6 +248,9 @@ def set_open_circuit_potential_submodel(self):
                     ocp_model = ocp_submodels.SingleOpenCircuitPotential
                 elif ocp_option == "current sigmoid":
                     ocp_model = ocp_submodels.CurrentSigmoidOpenCircuitPotential
+                elif ocp_option == "Wycisk":
+                    pybamm.citations.register("Wycisk2022")
+                    ocp_model = ocp_submodels.WyciskOpenCircuitPotential
                 elif ocp_option == "MSMR":
                     ocp_model = ocp_submodels.MSMROpenCircuitPotential
                 self.submodels[f"{domain} {phase} open-circuit potential"] = ocp_model(
diff --git a/pybamm/models/submodels/interface/open_circuit_potential/__init__.py b/pybamm/models/submodels/interface/open_circuit_potential/__init__.py
index 5f8a409bba..618a130f65 100644
--- a/pybamm/models/submodels/interface/open_circuit_potential/__init__.py
+++ b/pybamm/models/submodels/interface/open_circuit_potential/__init__.py
@@ -2,3 +2,4 @@
 from .single_ocp import SingleOpenCircuitPotential
 from .current_sigmoid_ocp import CurrentSigmoidOpenCircuitPotential
 from .msmr_ocp import MSMROpenCircuitPotential
+from .wycisk_ocp import WyciskOpenCircuitPotential
diff --git a/pybamm/models/submodels/interface/open_circuit_potential/wycisk_ocp.py b/pybamm/models/submodels/interface/open_circuit_potential/wycisk_ocp.py
new file mode 100644
index 0000000000..2a79c78a48
--- /dev/null
+++ b/pybamm/models/submodels/interface/open_circuit_potential/wycisk_ocp.py
@@ -0,0 +1,176 @@
+#
+# from Wycisk 2022
+#
+import pybamm
+from . import BaseOpenCircuitPotential
+
+
+class WyciskOpenCircuitPotential(BaseOpenCircuitPotential):
+    """
+    Class for open-circuit potential with hysteresis based on the approach outlined by Wycisk :footcite:t:'Wycisk2022'.
+    This approach employs a differential capacity hysteresis state variable. The decay and switching of the hysteresis state
+    is tunable via two additional parameters. The hysteresis state is updated based on the current and the differential capacity.
+    """
+
+    def get_fundamental_variables(self):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+        h = pybamm.Variable(
+            f"{Domain} electrode {phase_name}hysteresis state",
+            domains={
+                "primary": f"{domain} electrode",
+                "secondary": "current collector",
+            },
+        )
+        return {
+            f"{Domain} electrode {phase_name}hysteresis state": h,
+        }
+
+    def get_coupled_variables(self, variables):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+        phase = self.phase
+
+        if self.reaction == "lithium-ion main":
+            T = variables[f"{Domain} electrode temperature [K]"]
+            h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+
+            # For "particle-size distribution" models, take distribution version
+            # of c_s_surf that depends on particle size.
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "distribution":
+                sto_surf = variables[
+                    f"{Domain} {phase_name}particle surface stoichiometry distribution"
+                ]
+                # If variable was broadcast, take only the orphan
+                if isinstance(sto_surf, pybamm.Broadcast) and isinstance(
+                    T, pybamm.Broadcast
+                ):
+                    sto_surf = sto_surf.orphans[0]
+                    T = T.orphans[0]
+                T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"])
+                h = pybamm.PrimaryBroadcast(h, [f"{domain} {phase_name}particle size"])
+            else:
+                sto_surf = variables[
+                    f"{Domain} {phase_name}particle surface stoichiometry"
+                ]
+                # If variable was broadcast, take only the orphan
+                if isinstance(sto_surf, pybamm.Broadcast) and isinstance(
+                    T, pybamm.Broadcast
+                ):
+                    sto_surf = sto_surf.orphans[0]
+                    T = T.orphans[0]
+
+            variables[
+                f"{Domain} electrode {phase_name}hysteresis state distribution"
+            ] = h
+
+            # Bulk OCP is from the average SOC and temperature
+            sto_bulk = variables[f"{Domain} electrode {phase_name}stoichiometry"]
+            c_scale = self.phase_param.c_max
+            variables[f"Total lithium in {phase} phase in {domain} electrode [mol]"] = (
+                sto_bulk * c_scale
+            )  # c_s_vol * L * A
+
+            ocp_surf_eq = self.phase_param.U(sto_surf, T)
+            variables[f"{Domain} electrode {phase_name}equilibrium OCP [V]"] = (
+                ocp_surf_eq
+            )
+
+            T_bulk = pybamm.xyz_average(pybamm.size_average(T))
+            ocp_bulk_eq = self.phase_param.U(sto_bulk, T_bulk)
+            variables[f"{Domain} electrode {phase_name}bulk equilibrium OCP [V]"] = (
+                ocp_bulk_eq
+            )
+
+            inputs = {f"{Domain} {phase_name}particle stoichiometry": sto_surf}
+            lith_ref = pybamm.FunctionParameter(
+                f"{self.phase_param.phase_prefactor}{Domain} electrode lithiation OCP [V]",
+                inputs,
+            )
+            delith_ref = pybamm.FunctionParameter(
+                f"{self.phase_param.phase_prefactor}{Domain} electrode delithiation OCP [V]",
+                inputs,
+            )
+            H = abs(delith_ref - lith_ref) / 2
+            variables[f"{Domain} electrode {phase_name}OCP hysteresis [V]"] = H
+
+            # determine dQ/dU
+            if phase_name == "":
+                Q_mag = variables[f"{Domain} electrode capacity [A.h]"]
+            else:
+                Q_mag = variables[
+                    f"{Domain} electrode {phase_name}phase capacity [A.h]"
+                ]
+
+            dU = self.phase_param.U(sto_surf, T_bulk).diff(sto_surf)
+            dQdU = Q_mag / dU
+            variables[
+                f"{Domain} electrode {phase_name}differential capacity [A.s.V-1]"
+            ] = dQdU
+
+            H_x_av = pybamm.x_average(H)
+            h_x_av = pybamm.x_average(h)
+            variables[f"X-averaged {domain} electrode {phase_name}hysteresis state"] = (
+                h_x_av
+            )
+
+            # check if psd
+            if domain_options["particle size"] == "distribution":
+                # should always be true
+                if f"{domain} particle size" in sto_surf.domains["primary"]:
+                    # check if MPM Model
+                    if "current collector" in sto_surf.domains["secondary"]:
+                        ocp_surf = ocp_surf_eq + H_x_av * h_x_av
+                    # must be DFN with PSD model
+                    elif (
+                        f"{domain} electrode" in sto_surf.domains["secondary"]
+                        or f"{domain} {phase_name}particle size"
+                        in sto_surf.domains["primary"]
+                    ):
+                        ocp_surf = ocp_surf_eq + H * h
+            # must not be a psd
+            else:
+                ocp_surf = ocp_surf_eq + H * h
+
+            H_s_av = pybamm.size_average(H_x_av)
+            h_s_av = pybamm.size_average(h_x_av)
+
+            ocp_bulk = ocp_bulk_eq + H_s_av * h_s_av
+
+            dUdT = self.phase_param.dUdT(sto_surf)
+
+        variables.update(self._get_standard_ocp_variables(ocp_surf, ocp_bulk, dUdT))
+        return variables
+
+    def set_rhs(self, variables):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+
+        current = variables[
+            f"{Domain} electrode {phase_name}interfacial current density [A.m-2]"
+        ]
+        # check if composite or not
+        if phase_name != "":
+            Q_cell = variables[f"{Domain} electrode {phase_name}phase capacity [A.h]"]
+        else:
+            Q_cell = variables[f"{Domain} electrode capacity [A.h]"]
+
+        dQdU = variables[
+            f"{Domain} electrode {phase_name}differential capacity [A.s.V-1]"
+        ]
+        dQdU = dQdU.orphans[0]
+        K = self.phase_param.hysteresis_decay
+        K_x = self.phase_param.hysteresis_switch
+        h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+
+        dhdt = (
+            K * (current / (Q_cell * (dQdU**K_x))) * (1 - pybamm.sign(current) * h)
+        )  #! current is backwards for a halfcell
+        self.rhs[h] = dhdt
+
+    def set_initial_conditions(self, variables):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+        h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+        self.initial_conditions[h] = pybamm.Scalar(0)
diff --git a/pybamm/parameters/lithium_ion_parameters.py b/pybamm/parameters/lithium_ion_parameters.py
index 7ef032f7f1..72eb6534da 100644
--- a/pybamm/parameters/lithium_ion_parameters.py
+++ b/pybamm/parameters/lithium_ion_parameters.py
@@ -512,6 +512,14 @@ def _set_parameters(self):
         eps_c_init_av = pybamm.xyz_average(
             self.epsilon_s * pybamm.r_average(self.c_init)
         )
+        # if self.options['open-circuit potential'] == 'Plett':
+        self.hysteresis_decay = pybamm.Parameter(
+            f"{pref}{Domain} particle hysteresis decay rate"
+        )
+        self.hysteresis_switch = pybamm.Parameter(
+            f"{pref}{Domain} particle hysteresis switching factor"
+        )
+        self.h_init = pybamm.Scalar(0)
 
         if self.options["open-circuit potential"] != "MSMR":
             self.U_init = self.U(self.sto_init_av, main.T_init)
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
index e217a11d75..82d228badb 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
@@ -68,6 +68,28 @@ def test_current_sigmoid_ocp(self):
         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
         modeltest.test_all(skip_output_tests=True)
 
+    def test_wycisk_ocp(self):
+        options = {"open-circuit potential": ("Wycisk", "single")}
+        model = pybamm.lithium_ion.MPM(options)
+        parameter_values = pybamm.ParameterValues("Chen2020")
+        parameter_values = pybamm.get_size_distribution_parameters(parameter_values)
+        parameter_values.update(
+            {
+                "Negative electrode lithiation OCP [V]"
+                "": lambda sto: parameter_values["Negative electrode OCP [V]"](sto)
+                - 0.1,
+                "Negative electrode delithiation OCP [V]"
+                "": lambda sto: parameter_values["Negative electrode OCP [V]"](sto)
+                + 0.1,
+                "Negative particle hysteresis decay rate": 1,
+                "Negative particle hysteresis switching factor": 1,
+                # "Negative electrode OCP hysteresis [V]": lambda sto: 1,
+            },
+            check_already_exists=False,
+        )
+        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+        modeltest.test_all(skip_output_tests=True)
+
     def test_voltage_control(self):
         options = {"operating mode": "voltage"}
         model = pybamm.lithium_ion.MPM(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index c56cd2304c..5ff5027c1a 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -33,7 +33,7 @@
 'lithium plating porosity change': 'false' (possible: ['false', 'true'])
 'loss of active material': 'stress-driven' (possible: ['none', 'stress-driven', 'reaction-driven', 'current-driven', 'stress and reaction-driven'])
 'number of MSMR reactions': 'none' (possible: ['none'])
-'open-circuit potential': 'single' (possible: ['single', 'current sigmoid', 'MSMR'])
+'open-circuit potential': 'single' (possible: ['single', 'current sigmoid', 'MSMR', 'Wycisk'])
 'operating mode': 'current' (possible: ['current', 'voltage', 'power', 'differential power', 'explicit power', 'resistance', 'differential resistance', 'explicit resistance', 'CCCV'])
 'particle': 'Fickian diffusion' (possible: ['Fickian diffusion', 'fast diffusion', 'uniform profile', 'quadratic profile', 'quartic profile', 'MSMR'])
 'particle mechanics': 'swelling only' (possible: ['none', 'swelling only', 'swelling and cracking'])
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
index de7d2bfe83..4302681dc2 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
@@ -406,6 +406,10 @@ def test_well_posed_current_sigmoid_ocp(self):
         options = {"open-circuit potential": "current sigmoid"}
         self.check_well_posedness(options)
 
+    def test_well_posed_wycisk_ocp(self):
+        options = {"open-circuit potential": "Wycisk"}
+        self.check_well_posedness(options)
+
     def test_well_posed_msmr(self):
         options = {
             "open-circuit potential": "MSMR",
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index d7e95247e0..20fc69e541 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -35,6 +35,20 @@ def test_well_posed_current_sigmoid_ocp_with_psd(self):
         }
         self.check_well_posedness(options)
 
+    def test_well_posed_wycisk_ocp_with_psd(self):
+        options = {
+            "open-circuit potential": "Wycisk",
+            "particle size": "distribution",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_wycisk_ocp_with_composite(self):
+        options = {
+            "open-circuit potential": (("Wycisk", "single"), "single"),
+            "particle phases": ("2", "1"),
+        }
+        self.check_well_posedness(options)
+
     def test_well_posed_external_circuit_explicit_power(self):
         options = {"operating mode": "explicit power"}
         self.check_well_posedness(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
index 88049c0c63..389aa55849 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py
@@ -116,6 +116,13 @@ def test_msmr(self):
         model = pybamm.lithium_ion.MPM(options)
         model.check_well_posedness()
 
+    def test_wycisk_ocp(self):
+        options = {
+            "open-circuit potential": "Wycisk",
+        }
+        model = pybamm.lithium_ion.MPM(options)
+        model.check_well_posedness()
+
 
 class TestMPMExternalCircuits(TestCase):
     def test_well_posed_voltage(self):