diff --git a/content/02_basic_simpy.ipynb b/content/02_basic_simpy.ipynb
index b2887a8..488db42 100644
--- a/content/02_basic_simpy.ipynb
+++ b/content/02_basic_simpy.ipynb
@@ -230,7 +230,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.9"
+   "version": "3.11.11"
   }
  },
  "nbformat": 4,
diff --git a/content/11_oop_111.ipynb b/content/11_oop_111.ipynb
new file mode 100644
index 0000000..7ede9c8
--- /dev/null
+++ b/content/11_oop_111.ipynb
@@ -0,0 +1,770 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "de4dd3cf-09d1-45b9-b8fc-784edf3b0819",
+   "metadata": {},
+   "source": [
+    "# The full urgent care centre call model\n",
+    "\n",
+    "**We can now update the model to include the final call centre logic!**\n",
+    "\n",
+    "After a patient has spoken to a call operator their priority is triaged.  It is estimated that 40% of patients require a callback from a nurse.  There are 10 nurses available.  A nurse patient consultation has a Uniform distribution lasting between 10 and 20 minutes.\n",
+    "\n",
+    "> ⏰ Some call centres run 24/7 while others are open for specified time window during the day.  So we need to cope with **both terminating and non-terminating systems**.\n",
+    "\n",
+    "![model image](img/full_model.png \"Urgent care call centre\")\n",
+    "\n",
+    "**Modifications needed**\n",
+    "\n",
+    "* Add new default variables for the nurse consultation parameters\n",
+    "* Add new decision variables to `Experiment` for the no. nurses and the consultation distribution.\n",
+    "* Create a second `simpy.Resource` called `nurses` and add it to the simulation model.\n",
+    "* Create the nurse consultation process\n",
+    "* Modify the logic of `service` so that % of patients are called back.\n",
+    "* Collect results and estimate the waiting time for a nurse consultation and nurse utilisation.\n",
+    "* **Bonus:** Add an *optional* warm-up period event."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "57369b8c-451d-4d27-8ee0-c98d44d4e4e2",
+   "metadata": {},
+   "source": [
+    "## 1. Standard Imports"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "09f761a3-c100-4dc8-a3c9-30a5cd101c0d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import simpy\n",
+    "import itertools\n",
+    "\n",
+    "from abc import ABC, abstractmethod\n",
+    "from typing import Optional, Generator\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "d038d39b-84a3-4294-9e74-eae2192c437f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## 2. Model imports\n",
+    "\n",
+    "from model_utility import (\n",
+    "    Distribution,\n",
+    "    Exponential,\n",
+    "    Bernoulli,\n",
+    "    Uniform,\n",
+    "    Triangular,\n",
+    "    Experiment,\n",
+    "    WARM_UP_PERIOD,\n",
+    "    RESULTS_COLLECTION_PERIOD,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "08f6e5b4-ecf0-4da4-8c98-1b58f74e506e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def trace(msg):\n",
+    "    \"\"\"\n",
+    "    Turing printing of events on and off.\n",
+    "\n",
+    "    Params:\n",
+    "    -------\n",
+    "    msg: str\n",
+    "        string to print to screen.\n",
+    "    \"\"\"\n",
+    "    if TRACE:\n",
+    "        print(msg)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "33f7e4b8-46d1-4566-8307-2ad945455225",
+   "metadata": {},
+   "source": [
+    "## 3. Modified model code\n",
+    "\n",
+    "We will modify the model code and logic that we have already developed to include a nurse consultation for a proportion of the callers.  We create a new function called `nurse_consultation` that contains all the logic. We also need to modify the `service` function so that a proportion of calls are sent to the nurse consultation process.  "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "d4c51bf1-54a0-498e-97db-03273b7cf085",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class SimplePathwayDelay(ABC):\n",
+    "    \"\"\"\n",
+    "    Abstract base class representing a simple resource bound delay in a pathway \n",
+    "    simulation using SimPy.\n",
+    "    \n",
+    "    This class encapsulates the logic for a service activity within a pathway,\n",
+    "    including resource request, service duration, and custom logic before \n",
+    "    and after service.\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    def __init__(self, \n",
+    "                 identifier: int,                 \n",
+    "                 env: simpy.Environment, \n",
+    "                 resource: simpy.Resource,\n",
+    "                 delay_dist: Distribution,\n",
+    "                 args: Experiment,\n",
+    "                 trace_dp: Optional[int] = 2\n",
+    "    ) -> None:\n",
+    "        \"\"\"\n",
+    "        Initialize the pathway delay.\n",
+    "        \n",
+    "        Parameters:\n",
+    "        ----------            \n",
+    "        identifier: int\n",
+    "            Unique identifier for this delay\n",
+    "            \n",
+    "        env: simpy.Environment\n",
+    "            The simulation environment in which the delay process\n",
+    "            is running.\n",
+    "            \n",
+    "        resource: simpy.resource\n",
+    "            The simpy resource to be requested\n",
+    "\n",
+    "        delay_dist: Distribution \n",
+    "            Used for sampling the process activity duration. Protocol pattern\n",
+    "            must implement the sample method\n",
+    "            \n",
+    "        args: Experiment \n",
+    "            All model parameters\n",
+    "\n",
+    "        trace_dp: int, optional (default=2)\n",
+    "            The number of decimal places that the time is displayed to when\n",
+    "            sending a message to trace.\n",
+    "        \"\"\"\n",
+    "        self.identifier = identifier\n",
+    "        self.env = env\n",
+    "        self.resource = resource\n",
+    "        self.delay_dist = delay_dist\n",
+    "        self.args = args\n",
+    "        self.trace_dp = trace_dp\n",
+    "        self.waiting_time = -np.inf\n",
+    "        self.service_time = -np.inf\n",
+    "        self.total_time = -np.inf\n",
+    "            \n",
+    "    @abstractmethod\n",
+    "    def before_service_logic(self) -> None:\n",
+    "        \"\"\"Simulation logic prior to the process entering a service activity\"\"\"\n",
+    "        pass\n",
+    "\n",
+    "    @abstractmethod\n",
+    "    def after_service_logic(self) -> None:\n",
+    "        \"\"\"Simulation logic after a service activity has been completed\"\"\"\n",
+    "        pass\n",
+    "\n",
+    "    def service_logic(self) -> Generator:\n",
+    "        \"\"\"\n",
+    "        Simulate the service process, including resource request and activity duration.\n",
+    "        \n",
+    "        Returns:\n",
+    "        --------\n",
+    "            A generator that can be yielded in a SimPy environment\n",
+    "        \"\"\"\n",
+    "        # Record the start time for waiting\n",
+    "        start_wait = self.env.now\n",
+    "        \n",
+    "        # request a resource\n",
+    "        with self.resource.request() as req:\n",
+    "            yield req\n",
+    "    \n",
+    "            # record the waiting time for resource\n",
+    "            self.waiting_time = self.env.now - start_wait\n",
+    "                \n",
+    "            # sample the duration of the activity\n",
+    "            self.service_time = self.delay_dist.sample()\n",
+    "            \n",
+    "            # schedule process to begin again after activity duration\n",
+    "            yield self.env.timeout(self.service_time)\n",
+    "    \n",
+    "    def process(self) -> Generator:\n",
+    "        \"\"\"\n",
+    "        Run the complete service process, including before and after hooks.\n",
+    "        \n",
+    "        Returns:\n",
+    "            A generator that can be yielded in a SimPy environment\n",
+    "        \"\"\"\n",
+    "        start_time = self.env.now\n",
+    "        self.before_service_logic()\n",
+    "        yield from self.service_logic()\n",
+    "        self.after_service_logic()\n",
+    "        self.total_time = self.env.now - start_time\n",
+    "        \n",
+    "    def get_process(self) -> simpy.Process:\n",
+    "        \"\"\"\n",
+    "        Get a SimPy process for this delay that can be started in the \n",
+    "        environment.\n",
+    "        \n",
+    "        Returns:\n",
+    "        -------\n",
+    "            A SimPy process\n",
+    "        \"\"\"\n",
+    "        return self.env.process(self.process())\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "292f062b-676f-45a7-a744-163aa22011a8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class NurseConsultationPathway(SimplePathwayDelay):\n",
+    "    \"\"\"\n",
+    "    Simulates the wait for a consultation with a nurse on the phone 🩺\n",
+    "\n",
+    "    Inherits from SimplyPathwayDelay and hence follows this pattern:\n",
+    "    \n",
+    "    1. request and wait for a nurse resource\n",
+    "    2. phone consultation delay \n",
+    "    3. release nurse and exit system\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    def __init__(\n",
+    "        self, \n",
+    "        identifier: int, \n",
+    "        env: simpy.Environment, \n",
+    "        args: Experiment\n",
+    "    ) -> None:\n",
+    "        \"\"\"\n",
+    "        NurseConsultationPathway\n",
+    "\n",
+    "        Parameters:\n",
+    "        -----------        \n",
+    "        identifier: int\n",
+    "            Unique identifier for this delay\n",
+    "            \n",
+    "        env: simpy.Environment\n",
+    "            The simulation environment in which the delay process\n",
+    "            is running.\n",
+    "            \n",
+    "        args: Experiment \n",
+    "            All model parameters\n",
+    "        \"\"\"\n",
+    "        super().__init__(\n",
+    "            identifier=identifier,                 \n",
+    "            env=env,\n",
+    "            resource=args.nurses,\n",
+    "            delay_dist=args.nurse_dist,\n",
+    "            args=args)\n",
+    "        \n",
+    "        self.activity_name = \"nurse consultation\"\n",
+    "\n",
+    "    def before_service(self) -> None:\n",
+    "        \"\"\"Simulation logic prior to the process entering a service activity\"\"\"\n",
+    "        trace(\n",
+    "            f\"{self.env.now:.{self.trace_dp}f}: Patient {self.identifier}\"\n",
+    "            + f\"waiting for nurse call back\"\n",
+    "        )\n",
+    "\n",
+    "    def after_service(self) -> None:\n",
+    "        \"\"\"Simulation logic after service has completed.\"\"\"\n",
+    "        # log metrics for process\n",
+    "        self.args.results[\"nurse_waiting_times\"].append(self.waiting_time)\n",
+    "        self.args.results[\"total_nurse_call_duration\"] += self.service_time\n",
+    "\n",
+    "        trace(\n",
+    "            f\"{self.env.now:.{self.trace_dp}f}: {self.activity_name}\"\n",
+    "            + f\" for {self.identifier} completed.\"\n",
+    "        )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "56fa8b72-c511-44a8-99d4-d958ef256821",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class CallOperatorProcess(SimplePathwayDelay):\n",
+    "    \"\"\"\n",
+    "    Encapsulates a Call Operator Process 💬\n",
+    "    \n",
+    "    This is a SimplePathwayDelay i.e. a FIFO resource bound delay.\n",
+    "\n",
+    "    1. Request a call operator\n",
+    "    2. Simulate a delay while call operator triage takes place\n",
+    "    3. Release the resource\n",
+    "    4. Simulate the need for a nurse call back\n",
+    "    5. If needed launch an instance of NurseConsultationPathway\n",
+    "    \"\"\"\n",
+    "    def __init__(\n",
+    "        self, \n",
+    "        identifier: int, \n",
+    "        env: simpy.Environment, \n",
+    "        args: Experiment\n",
+    "    ) -> None:\n",
+    "        \"\"\"\n",
+    "        CallOperatorProcess\n",
+    "\n",
+    "        Parameters:\n",
+    "        -----------        \n",
+    "        identifier: int\n",
+    "            Unique identifier for this delay\n",
+    "            \n",
+    "        env: simpy.Environment\n",
+    "            The simulation environment in which the delay process\n",
+    "            is running.\n",
+    "            \n",
+    "        args: Experiment \n",
+    "            All model parameters\n",
+    "        \"\"\"\n",
+    "        \n",
+    "        super().__init__(\n",
+    "            identifier=identifier,                 \n",
+    "            env=env,\n",
+    "            resource=args.operators,\n",
+    "            delay_dist=args.call_dist,\n",
+    "            args=args)\n",
+    "        \n",
+    "        self.activity_name = \"Call operator service\"\n",
+    "\n",
+    "    def before_service(self) -> None:\n",
+    "        \"\"\"Simulation logic prior to the process entering a service activity\"\"\"\n",
+    "        trace(\n",
+    "            f\"{self.env.now:.{self.trace_dp}f}: Patient {self.identifier}\"\n",
+    "            + f\" waiting for nurse call back.\"\n",
+    "        )\n",
+    "\n",
+    "    def after_service(self) -> None:\n",
+    "        \"\"\"Simulation logic after a service has completed.\n",
+    "        \n",
+    "        1. record process metrics\n",
+    "        2. simulate nurse call back process\n",
+    "        \n",
+    "        \"\"\"\n",
+    "        # log metrics for process\n",
+    "        self.args.results[\"waiting_times\"].append(self.waiting_time)\n",
+    "        self.args.results[\"total_call_duration\"] += self.service_time\n",
+    "\n",
+    "        trace(\n",
+    "            f\"{self.env.now:.{self.trace_dp}f}: {self.activity_name}\"\n",
+    "            + f\" for {self.identifier} completed.\"\n",
+    "        )\n",
+    "\n",
+    "        # nurse consultation?\n",
+    "        self.next_process()\n",
+    "\n",
+    "    def next_process(self) -> None:\n",
+    "        \"\"\"Simulate if a nurse call back is required.\"\"\"\n",
+    "        \n",
+    "        callback_patient = self.args.callback_dist.sample()\n",
+    "        \n",
+    "        if callback_patient:\n",
+    "            nurse_pathway = NurseConsultationPathway(\n",
+    "                self.identifier, \n",
+    "                self.env, \n",
+    "                self.args\n",
+    "            )\n",
+    "            self.env.process(nurse_pathway.process())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "838b63ac-8add-442a-b25b-ed9c7e5fdfd9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class UrgentCallCentreModel:\n",
+    "    \"\"\"\n",
+    "    Urgent care call centre model\n",
+    "    \"\"\"\n",
+    "    def __init__(\n",
+    "        self,\n",
+    "        env: simpy.Environment, \n",
+    "        args: Experiment\n",
+    "    ) -> None:\n",
+    "        \"\"\"\n",
+    "        UrgentCallCentreModel\n",
+    "\n",
+    "        Parameters:\n",
+    "        ------\n",
+    "        env: simpy.Environment\n",
+    "            The simpy environment for the simulation\n",
+    "    \n",
+    "        args: Experiment\n",
+    "            The settings and input parameters for the simulation.\n",
+    "        \n",
+    "        \"\"\"\n",
+    "        self.env = env\n",
+    "        self.args = args\n",
+    "\n",
+    "        self._init_resources()\n",
+    "\n",
+    "    def _init_resources(self) -> None:\n",
+    "        \"\"\"\n",
+    "        Init the number of resources\n",
+    "        and store in the arguments container object\n",
+    "        \n",
+    "        Resource list:\n",
+    "        1. call operators\n",
+    "        2. nurses\n",
+    "        \"\"\"\n",
+    "        # create the environment object.\n",
+    "        self.args.operators = simpy.Resource(self.env, self.args.n_operators)\n",
+    "    \n",
+    "        # create the nurses resource\n",
+    "        self.args.nurses = simpy.Resource(self.env, self.args.n_nurses)\n",
+    "\n",
+    "    def arrivals_generator(self) -> Generator:\n",
+    "        \"\"\"\n",
+    "        📞 Call arrivals generator\n",
+    "        \"\"\"\n",
+    "        # use itertools as it provides an infinite loop\n",
+    "        # with a counter variable that we can use for unique Ids\n",
+    "        for caller_count in itertools.count(start=1):\n",
+    "    \n",
+    "            # rhe sample distribution is defined by the experiment.\n",
+    "            inter_arrival_time = self.args.arrival_dist.sample()\n",
+    "            yield self.env.timeout(inter_arrival_time)\n",
+    "    \n",
+    "            trace(f\"{self.env.now:.2f}: call arrives.\")\n",
+    "    \n",
+    "            # create a service process\n",
+    "            operator_pathway = CallOperatorProcess(\n",
+    "                caller_count, \n",
+    "                self.env, \n",
+    "                self.args\n",
+    "            )\n",
+    "\n",
+    "            # initialise pathway as simpy process\n",
+    "            self.env.process(operator_pathway.process())\n",
+    "\n",
+    "    def warmup_complete(self, warm_up_period: float) -> Generator:\n",
+    "        \"\"\"\n",
+    "        🥵 End of warm-up period Generator. \n",
+    "        Generates a single event that occurs at the end of the warm-up period\n",
+    "        Used to reset results collection variables.\n",
+    "    \n",
+    "        Parameters:\n",
+    "        ----------\n",
+    "        warm_up_period: float\n",
+    "            Duration of warm-up period in simultion time units\n",
+    "        \"\"\"\n",
+    "        yield self.env.timeout(warm_up_period)\n",
+    "        trace(f\"{self.env.now:.2f}: Warm up complete.\")\n",
+    "\n",
+    "        # reset KPIs\n",
+    "        self.args.init_results_variables()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "161ce6c8-b50b-457c-8f9e-4837edb07933",
+   "metadata": {},
+   "source": [
+    "## 5. Model wrapper functions\n",
+    "\n",
+    "Modifications to make to the `single_run` function:\n",
+    "\n",
+    "1. Add a warm-up parameters called `wu_period`\n",
+    "1. Create and the nurses resource to the experiment\n",
+    "2. Schedule the `warm_up_complete` process.\n",
+    "3. After the simulation is complete calculate the mean waiting time and mean utilisation for nurses."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "1bde92ab-c010-4fa9-90cf-ebb715537ae9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def single_run(\n",
+    "    experiment: Experiment, \n",
+    "    rep: Optional[int] = 0,\n",
+    "    wu_period: Optional[float] = WARM_UP_PERIOD, \n",
+    "    rc_period: Optional[float] = RESULTS_COLLECTION_PERIOD\n",
+    ") -> dict:\n",
+    "    \"\"\"\n",
+    "    Perform a single run of the model and return the results\n",
+    "\n",
+    "    Parameters:\n",
+    "    -----------\n",
+    "\n",
+    "    experiment: Experiment\n",
+    "        The experiment/paramaters to use with model\n",
+    "\n",
+    "    rep: int\n",
+    "        The replication number.\n",
+    "\n",
+    "    wu_period: float, optional (default=WARM_UP_PERIOD)\n",
+    "        The initial transient period of the simulation\n",
+    "        Results from this period are removed from final computations.\n",
+    "\n",
+    "    rc_period: float, optional (default=RESULTS_COLLECTION_PERIOD)\n",
+    "        The run length of the model following warm up where results are\n",
+    "        collected.\n",
+    "    \"\"\"\n",
+    "\n",
+    "    # results dictionary.  Each KPI is a new entry.\n",
+    "    run_results = {}\n",
+    "\n",
+    "    # reset all results variables to zero and empty\n",
+    "    experiment.init_results_variables()\n",
+    "\n",
+    "    # set random number set to the replication no.\n",
+    "    # this controls sampling for the run.\n",
+    "    experiment.set_random_no_set(rep)\n",
+    "\n",
+    "    # environment is (re)created inside single run\n",
+    "    env = simpy.Environment()\n",
+    "    model = UrgentCallCentreModel(env, experiment)\n",
+    "    \n",
+    "    # we pass the experiment to the arrivals generator\n",
+    "    env.process(model.arrivals_generator())\n",
+    "\n",
+    "    # add warm-up period event\n",
+    "    env.process(model.warmup_complete(wu_period))\n",
+    "\n",
+    "    # run for warm-up + results collection period\n",
+    "    env.run(until=wu_period + rc_period)\n",
+    "\n",
+    "    # end of run results: calculate mean waiting time\n",
+    "    run_results[\"01_mean_waiting_time\"] = np.mean(\n",
+    "        experiment.results[\"waiting_times\"]\n",
+    "    )\n",
+    "\n",
+    "    # end of run results: calculate mean operator utilisation\n",
+    "    run_results[\"02_operator_util\"] = (\n",
+    "        experiment.results[\"total_call_duration\"]\n",
+    "        / (rc_period * experiment.n_operators)\n",
+    "    ) * 100.0\n",
+    "\n",
+    "    # end of run results: nurse waiting time\n",
+    "    run_results[\"03_mean_nurse_waiting_time\"] = np.mean(\n",
+    "        experiment.results[\"nurse_waiting_times\"]\n",
+    "    )\n",
+    "\n",
+    "    # end of run results: calculate mean nurse utilisation\n",
+    "    run_results[\"04_nurse_util\"] = (\n",
+    "        experiment.results[\"total_nurse_call_duration\"]\n",
+    "        / (rc_period * experiment.n_nurses)\n",
+    "    ) * 100.0\n",
+    "\n",
+    "    # return the results from the run of the model\n",
+    "    return run_results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "c6997bce-ae2e-449a-bf9f-6ba09acb6d55",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def multiple_replications(\n",
+    "    experiment: Experiment,\n",
+    "    wu_period: Optional[float] = WARM_UP_PERIOD, \n",
+    "    rc_period: Optional[float] = RESULTS_COLLECTION_PERIOD,\n",
+    "    n_reps: Optional[int] = 5,\n",
+    ") -> pd.DataFrame:\n",
+    "    \"\"\"\n",
+    "    Perform multiple replications of the model.\n",
+    "\n",
+    "    Params:\n",
+    "    ------\n",
+    "    experiment: Experiment\n",
+    "        The experiment/paramaters to use with model\n",
+    "\n",
+    "    rc_period: float, optional (default=DEFAULT_RESULTS_COLLECTION_PERIOD)\n",
+    "        results collection period.\n",
+    "        the number of minutes to run the model to collect results\n",
+    "\n",
+    "    n_reps: int, optional (default=5)\n",
+    "        Number of independent replications to run.\n",
+    "\n",
+    "    Returns:\n",
+    "    --------\n",
+    "    pandas.DataFrame\n",
+    "    \"\"\"\n",
+    "\n",
+    "    # loop over single run to generate results dicts in a python list.\n",
+    "    results = [\n",
+    "        single_run(experiment, rep, wu_period, rc_period) \n",
+    "        for rep in range(n_reps)\n",
+    "    ]\n",
+    "\n",
+    "    # format and return results in a dataframe\n",
+    "    df_results = pd.DataFrame(results)\n",
+    "    df_results.index = np.arange(1, len(df_results) + 1)\n",
+    "    df_results.index.name = \"rep\"\n",
+    "    return df_results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "50fcc6ef-9e7d-4caf-8f62-fbe4e4e1d061",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>count</th>\n",
+       "      <th>mean</th>\n",
+       "      <th>std</th>\n",
+       "      <th>min</th>\n",
+       "      <th>25%</th>\n",
+       "      <th>50%</th>\n",
+       "      <th>75%</th>\n",
+       "      <th>max</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>01_mean_waiting_time</th>\n",
+       "      <td>5.0</td>\n",
+       "      <td>3.2</td>\n",
+       "      <td>0.9</td>\n",
+       "      <td>1.8</td>\n",
+       "      <td>3.0</td>\n",
+       "      <td>3.2</td>\n",
+       "      <td>3.9</td>\n",
+       "      <td>4.0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>02_operator_util</th>\n",
+       "      <td>5.0</td>\n",
+       "      <td>94.2</td>\n",
+       "      <td>1.3</td>\n",
+       "      <td>92.1</td>\n",
+       "      <td>94.1</td>\n",
+       "      <td>94.4</td>\n",
+       "      <td>94.5</td>\n",
+       "      <td>95.8</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>03_mean_nurse_waiting_time</th>\n",
+       "      <td>5.0</td>\n",
+       "      <td>3.0</td>\n",
+       "      <td>1.7</td>\n",
+       "      <td>1.6</td>\n",
+       "      <td>1.6</td>\n",
+       "      <td>2.0</td>\n",
+       "      <td>4.7</td>\n",
+       "      <td>4.9</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>04_nurse_util</th>\n",
+       "      <td>5.0</td>\n",
+       "      <td>88.3</td>\n",
+       "      <td>3.6</td>\n",
+       "      <td>83.7</td>\n",
+       "      <td>86.9</td>\n",
+       "      <td>87.9</td>\n",
+       "      <td>89.5</td>\n",
+       "      <td>93.6</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "                            count  mean  std   min   25%   50%   75%   max\n",
+       "01_mean_waiting_time          5.0   3.2  0.9   1.8   3.0   3.2   3.9   4.0\n",
+       "02_operator_util              5.0  94.2  1.3  92.1  94.1  94.4  94.5  95.8\n",
+       "03_mean_nurse_waiting_time    5.0   3.0  1.7   1.6   1.6   2.0   4.7   4.9\n",
+       "04_nurse_util                 5.0  88.3  3.6  83.7  86.9  87.9  89.5  93.6"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "TRACE = False\n",
+    "scenario = Experiment(n_nurses=15, nurse_call_high=30.0)\n",
+    "results = multiple_replications(scenario, wu_period=50.0)\n",
+    "results.describe().round(1).T"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "e03a228e-b642-405e-a9b2-ba383ee1fd8e",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'01_mean_waiting_time': 3.8923863988435095,\n",
+       " '02_operator_util': 94.05731922168263,\n",
+       " '03_mean_nurse_waiting_time': 4.945822358637601,\n",
+       " '04_nurse_util': 93.57384685068052}"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "TRACE = False\n",
+    "scenario = Experiment(n_nurses=15, nurse_call_high=30.0)\n",
+    "results = single_run(scenario, wu_period=50.0)\n",
+    "results"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/content/model_utility.py b/content/model_utility.py
new file mode 100644
index 0000000..9165fec
--- /dev/null
+++ b/content/model_utility.py
@@ -0,0 +1,359 @@
+import numpy as np
+import pandas as pd
+from typing import Protocol, Optional, Union, runtime_checkable
+from numpy.random import SeedSequence
+
+# Constants
+# default resources
+N_OPERATORS = 13
+N_NURSES = 10
+
+# default mean inter-arrival time (exp)
+MEAN_IAT = 60 / 100
+
+## default service time parameters (triangular)
+CALL_LOW = 5.0
+CALL_MODE = 7.0
+CALL_HIGH = 10.0
+
+# nurse uniform distribution parameters
+NURSE_CALL_LOW = 10.0
+NURSE_CALL_HIGH = 20.0
+
+# probability of a callback (parameter of Bernoulli)
+CHANCE_CALLBACK = 0.4
+
+# sampling settings - we now need 4 streams
+N_STREAMS = 4
+DEFAULT_RND_SET = 0
+
+# Boolean switch to simulation results as the model runs
+TRACE = False
+
+# run variables
+RESULTS_COLLECTION_PERIOD = 1000
+WARM_UP_PERIOD = 0
+
+# DISTRIBUTION CLASSES
+
+
+@runtime_checkable
+class Distribution(Protocol):
+    def sample(self, size: Optional[int] = None) -> Union[float, np.ndarray]:
+        """Generate random samples from this distribution.
+        
+        Parameters
+        ----------
+        size : Optional[int], default=None
+            If None, return a single sample.
+            If int, return that many samples in a 1D array.
+                 
+        Returns
+        -------
+        Union[float, np.ndarray]
+            A single sample (if size is None) or array of samples
+        """
+        pass
+
+
+class Bernoulli:
+    """
+    Convenience class for the Bernoulli distribution.
+    Packages up distribution parameters, seed and random generator.
+
+    The Bernoulli distribution is a special case of the binomial distribution
+    where a single trial is conducted.
+
+    Use the Bernoulli distribution to sample success or failure.
+    """
+
+    def __init__(self, p: float, random_seed: Optional[Union[int, SeedSequence]] = None):
+        """
+        Constructor
+
+        Parameters
+        ----------
+        p : float
+            Probability of drawing a 1
+            
+        random_seed : Optional[Union[int, SeedSequence]], default=None
+            A random seed to reproduce samples. If set to None then a unique
+            sample is created.
+        """
+        self.rand = np.random.default_rng(seed=random_seed)
+        self.p = p
+
+    def sample(self, size: Optional[int] = None) -> Union[float, np.ndarray]:
+        """
+        Generate a sample from the Bernoulli distribution.
+
+        Parameters
+        ----------
+        size : Optional[int], default=None
+            The number of samples to return. If size=None then a single
+            sample is returned.
+
+        Returns
+        -------
+        Union[float, np.ndarray]
+            Either a single sample (if size is None) or an array of samples
+        """
+        return self.rand.binomial(n=1, p=self.p, size=size)
+
+
+class Uniform:
+    """
+    Convenience class for the Uniform distribution.
+    Packages up distribution parameters, seed and random generator.
+    """
+
+    def __init__(self, low: float, high: float, random_seed: Optional[Union[int, SeedSequence]] = None):
+        """
+        Constructor
+
+        Parameters
+        ----------
+        low : float
+            Lower range of the uniform
+        
+        high : float
+            Upper range of the uniform
+        
+        random_seed : Optional[Union[int, SeedSequence]], default=None
+            A random seed to reproduce samples. If set to None then a unique
+            sample is created.
+        """
+        self.rand = np.random.default_rng(seed=random_seed)
+        self.low = low
+        self.high = high
+
+    def sample(self, size: Optional[int] = None) -> Union[float, np.ndarray]:
+        """
+        Generate a sample from the uniform distribution
+
+        Parameters
+        ----------
+        size : Optional[int], default=None
+            The number of samples to return. If size=None then a single
+            sample is returned.
+
+        Returns
+        -------
+        Union[float, np.ndarray]
+            Either a single sample (if size is None) or an array of samples
+        """
+        return self.rand.uniform(low=self.low, high=self.high, size=size)
+
+class Triangular:
+    """
+    Convenience class for the triangular distribution.
+    Packages up distribution parameters, seed and random generator.
+    """
+
+    def __init__(self, low: float, mode: float, high: float, 
+                 random_seed: Optional[Union[int, SeedSequence]] = None):
+        """
+        Constructor. Accepts and stores parameters of the triangular dist
+        and a random seed.
+
+        Parameters
+        ----------
+        low : float
+            The smallest values that can be sampled
+        
+        mode : float
+            The most frequently sampled value
+        
+        high : float
+            The highest value that can be sampled
+        
+        random_seed : Optional[Union[int, SeedSequence]], default=None
+            Used with params to create a series of repeatable samples.
+        """
+        self.rand = np.random.default_rng(seed=random_seed)
+        self.low = low
+        self.high = high
+        self.mode = mode
+
+    def sample(self, size: Optional[int] = None) -> Union[float, np.ndarray]:
+        """
+        Generate one or more samples from the triangular distribution
+
+        Parameters
+        ----------
+        size : Optional[int], default=None
+            The number of samples to return. If size=None then a single
+            sample is returned.
+
+        Returns
+        -------
+        Union[float, np.ndarray]
+            Either a single sample (if size is None) or an array of samples
+        """
+        return self.rand.triangular(self.low, self.mode, self.high, size=size)
+
+class Exponential:
+    """
+    Convenience class for the exponential distribution.
+    Packages up distribution parameters, seed and random generator.
+    """
+
+    def __init__(self, mean: float, random_seed: Optional[Union[int, SeedSequence]] = None):
+        """
+        Constructor
+
+        Parameters
+        ----------
+        mean : float
+            The mean of the exponential distribution
+        
+        random_seed : Optional[Union[int, SeedSequence]], default=None
+            A random seed to reproduce samples. If set to None then a unique
+            sample is created.
+        """
+        self.rand = np.random.default_rng(seed=random_seed)
+        self.mean = mean
+
+    def sample(self, size: Optional[int] = None) -> Union[float, np.ndarray]:
+        """
+        Generate a sample from the exponential distribution
+
+        Parameters
+        ----------
+        size : Optional[int], default=None
+            The number of samples to return. If size=None then a single
+            sample is returned.
+
+        Returns
+        -------
+        Union[float, np.ndarray]
+            Either a single sample (if size is None) or an array of samples
+        """
+        return self.rand.exponential(self.mean, size=size)
+
+
+# Experiment class
+
+class Experiment:
+    """
+    Encapsulates the concept of an experiment 🧪 with the urgent care
+    call centre simulation model.
+
+    An Experiment:
+    1. Contains a list of parameters that can be left as defaults or varied
+    2. Provides a place for the experimentor to record results of a run
+    3. Controls the set & streams of psuedo random numbers used in a run.
+
+    """
+
+    def __init__(
+        self,
+        random_number_set=DEFAULT_RND_SET,
+        n_streams=N_STREAMS,
+        n_operators=N_OPERATORS,
+        mean_iat=MEAN_IAT,
+        call_low=CALL_LOW,
+        call_mode=CALL_MODE,
+        call_high=CALL_HIGH,
+        n_nurses=N_NURSES,
+        chance_callback=CHANCE_CALLBACK,
+        nurse_call_low=NURSE_CALL_LOW,
+        nurse_call_high=NURSE_CALL_HIGH,
+    ):
+        """
+        The init method sets up our defaults.
+        """
+        # sampling
+        self.random_number_set = random_number_set
+        self.n_streams = n_streams
+
+        # store parameters for the run of the model
+        self.n_operators = n_operators
+        self.mean_iat = mean_iat
+        self.call_low = call_low
+        self.call_mode = call_mode
+        self.call_high = call_high
+
+        # resources: we must init resources after an Environment is created.
+        # But we will store a placeholder for transparency
+        self.operators = None
+
+        # nurse parameters
+        self.n_nurses = n_nurses
+        self.chance_callback = chance_callback
+        self.nurse_call_low = nurse_call_low
+        self.nurse_call_high = nurse_call_high
+
+        # nurse resources placeholder
+        self.nurses = None
+
+        # initialise results to zero
+        self.init_results_variables()
+
+        # initialise sampling objects
+        self.init_sampling()
+
+    def set_random_no_set(self, random_number_set):
+        """
+        Controls the random sampling
+        Parameters:
+        ----------
+        random_number_set: int
+            Used to control the set of pseudo random numbers used by
+            the distributions in the simulation.
+        """
+        self.random_number_set = random_number_set
+        self.init_sampling()
+
+    def init_sampling(self):
+        """
+        Create the distributions used by the model and initialise
+        the random seeds of each.
+        """
+        # produce n non-overlapping streams
+        seed_sequence = np.random.SeedSequence(self.random_number_set)
+        self.seeds = seed_sequence.spawn(self.n_streams)
+
+        # create distributions
+
+        # call inter-arrival times
+        self.arrival_dist = Exponential(
+            self.mean_iat, random_seed=self.seeds[0]
+        )
+
+        # duration of call triage
+        self.call_dist = Triangular(
+            self.call_low,
+            self.call_mode,
+            self.call_high,
+            random_seed=self.seeds[1],
+        )
+
+        # create the callback and nurse consultation distributions
+        self.callback_dist = Bernoulli(
+            self.chance_callback, random_seed=self.seeds[2]
+        )
+
+        self.nurse_dist = Uniform(
+            self.nurse_call_low,
+            self.nurse_call_high,
+            random_seed=self.seeds[3],
+        )
+
+    def init_results_variables(self):
+        """
+        Initialise all of the experiment variables used in results
+        collection.  This method is called at the start of each run
+        of the model
+        """
+        # variable used to store results of experiment
+        self.results = {}
+        self.results["waiting_times"] = []
+
+        # total operator usage time for utilisation calculation.
+        self.results["total_call_duration"] = 0.0
+
+        # nurse sub process results collection
+        self.results["nurse_waiting_times"] = []
+        self.results["total_nurse_call_duration"] = 0.0
+
