Hopfenberg

CHANGELOG.md

@@ -5,6 +5,8 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/)
 and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Added
+- Hopfenberg model
 ## [0.3] - 2025-12-08
 ### Added
 - Weibull model

README.md

@@ -150,6 +150,29 @@ where:
 3. Comparative studies
 4. Vivo predictions
 
+### Hopfenberg
+The Hopfenberg model describes drug release from surface-eroding polymers with a constant surface area. 
+It is particularly useful for modeling drug release from biodegradable polymers where the drug is uniformly distributed throughout the matrix.
+According to this model, the cumulative amount of drug released at time $t$ is given by:
+
+$$
+M_t = M_{\infty} \left(1 - \left(1 - \frac{k_0t}{c_0 a_0}\right)^n\right)
+$$
+
+where:
+- $M_t (mg)$ is the cumulative absolute amount of drug released at time $t
+- $M_{\infty} (mg)$ is the total amount of drug released at infinite time
+- $k_0 (\frac{mg}{mm^2 s})$ is the surface erosion rate constant
+- $c_0 (\frac{mg}{mm^3})$ is the initial drug concentration in the polymer matrix
+- $a_0 (mm)$ is the initial radius (for spheres) or half-thickness (for slabs) of the device
+- $n$ is the geometry-dependent exponent (1 for slabs, 2 for cylinders, 3 for spheres)
+
+#### Applications
+1. Biodegradable Implants
+2. Surface-eroding drug delivery systems
+3. Transdermal Patches
+4. Injectable depots
+
 ## Usage
 ### Zero-Order Model
 ```python
@@ -188,6 +211,17 @@ model.plot(show=True)
 ```
 <img src="https://github.com/openscilab/drux/raw/main/otherfiles/weibull_plot.png" alt="Weibull Plot">
 
+### Hopfenberg Model
+
+```python
+from drux import HopfenbergModel
+model = HopfenbergModel(M=1, k0=0.00067, c0=0.0374, a0=3.51, n=2)
+model.simulate(duration=100, time_step=1)
+model.plot(show=True)
+```
+<img src="https://github.com/openscilab/drux/raw/main/otherfiles/hopfenberg_plot.png" alt="Hopfenberg Plot">
+
+
 ## Issues & bug reports
 
 Just fill an issue and describe it. We'll check it ASAP! or send an email to [drux@openscilab.com](mailto:drux@openscilab.com "drux@openscilab.com"). 
@@ -209,6 +243,9 @@ You can also join our discord server
 <blockquote>5- S. Dash, "Kinetic modeling on drug release from controlled drug delivery systems," <i>Acta Poloniae Pharmaceutica</i>, 2010.</blockquote>
 <blockquote>6- K. H. Ramteke, P. A. Dighe, A. R. Kharat, S. V. Patil, <i>Mathematical models of drug dissolution: A review</i>, <i>Sch. Acad. J. Pharm.</i>, vol. 3, no. 5, pp. 388-396, 2014.</blockquote>
 <blockquote>7- C. Corsaro, G. Neri, A. M. Mezzasalma, and E. Fazio, "Weibull modeling of controlled drug release from Ag-PMA nanosystems," <i>Polymers</i>, vol. 13, no. 17, p. 2897, 2021.</blockquote>
+<blockquote>8- H. B. Hopfenberg, "Controlled release from erodible slabs, cylinders, and spheres," in <i>Controlled Release Polymeric Formulations</i>, ACS Symposium Series, vol. 33, pp. 26–32, 1976.</blockquote>
+<blockquote>9- H. V. Chavda, M. S. Patel, and C. N. Patel, "Preparation and in vitro evaluation of guar gum based triple-layer matrix tablet of diclofenac sodium," <i>Research in Pharmaceutical Sciences</i>, vol. 7, no. 1, pp. 57–64, Jan. 2012.</blockquote>
+
 
 ## Show your support
 ### Star this repo

drux/__init__.py

@@ -6,5 +6,6 @@
 from .zero_order import ZeroOrderModel, ZeroOrderParameters
 from .first_order import FirstOrderModel, FirstOrderParameters
 from .weibull import WeibullModel, WeibullParameters
+from .hopfenberg import HopfenbergModel, HopfenbergParameters
 
 __version__ = DRUX_VERSION

drux/hopfenberg.py

@@ -0,0 +1,92 @@
+# -*- coding: utf-8 -*-
+"""Drux Hopfenberg model implementation."""
+
+from .base_model import DrugReleaseModel
+from .messages import (
+    ERROR_INVALID_EROSION_CONSTANT,
+    ERROR_INVALID_INITIAL_RADIUS,
+    ERROR_INVALID_GEOMETRY_FACTOR,
+    ERROR_INVALID_CONCENTRATION,
+    ERROR_RELEASABLE_AMOUNT,
+)
+from dataclasses import dataclass
+
+
+@dataclass
+class HopfenbergParameters:
+    """
+    Parameters for the Hopfenberg model based on surface erosion.
+
+    Attributes:
+        k0 (float): Erosion rate constant (mg/(mm^2·s))
+        c0 (float): Initial drug concentration in the matrix (mg/mm^3)
+        a0 (float): Initial radius or half-thickness of the device (mm)
+        n (int): Geometry factor (1=slab, 2=cylinder, 3=sphere)
+    """
+
+    M: float
+    k0: float
+    c0: float
+    a0: float
+    n: int
+
+
+class HopfenbergModel(DrugReleaseModel):
+    """Simulator for the Hopfenberg drug release model for surface-eroding polymers."""
+
+    def __init__(self, M: float, k0: float, c0: float, a0: float, n: int) -> None:
+        """
+        Initialize the Hopfenberg model with the given parameters.
+
+        :param M: entire releasable amount of drug (normally M > 0) (mg)
+        :param k0: Erosion rate constant (mg/(mm^2·s))
+        :param c0: Initial drug concentration in the matrix (mg/mm^3)
+        :param a0: Initial radius or half-thickness of the device (mm)
+        :param n: Geometry factor (1=slab, 2=cylinder, 3=sphere)
+        """
+        super().__init__()
+        self._parameters = HopfenbergParameters(M=M, k0=k0, c0=c0, a0=a0, n=n)
+        self._plot_parameters["label"] = "Hopfenberg Model"
+
+    def __repr__(self):
+        """Return a string representation of the Hopfenberg model."""
+        return (
+            f"drux.HopfenbergModel(M={self._parameters.M}, k0={self._parameters.k0}, "
+            f"c0={self._parameters.c0}, a0={self._parameters.a0}, "
+            f"n={self._parameters.n})"
+        )
+
+    def _model_function(self, t: float) -> float:
+        """
+        Calculate the fractional drug release at time t using the Hopfenberg model.
+
+        Formula:
+        - Mt = M∞(1 - (1 - k0*t / (c0*a0))^n)
+
+        :param t: time (s)
+        :return: drug release
+        """
+        M = self._parameters.M
+        k0 = self._parameters.k0
+        c0 = self._parameters.c0
+        a0 = self._parameters.a0
+        n = self._parameters.n
+
+        inner_term = 1 - (k0 * t) / (c0 * a0)
+
+        Mt = M * (1 - (inner_term**n))
+
+        return Mt
+
+    def _validate_parameters(self) -> None:
+        """Validate the parameters of the Hopfenberg model."""
+        if self._parameters.M < 0:
+            raise ValueError(ERROR_RELEASABLE_AMOUNT)
+        if self._parameters.k0 < 0:
+            raise ValueError(ERROR_INVALID_EROSION_CONSTANT)
+        if self._parameters.c0 <= 0:
+            raise ValueError(ERROR_INVALID_CONCENTRATION)
+        if self._parameters.a0 <= 0:
+            raise ValueError(ERROR_INVALID_INITIAL_RADIUS)
+        if self._parameters.n not in (1, 2, 3):
+            raise ValueError(ERROR_INVALID_GEOMETRY_FACTOR)

drux/messages.py

@@ -36,6 +36,11 @@
 # Error messages for Weibull
 ERROR_WEIBULL_SCALE_PARAMETER = "Scale parameter (a) must be positive."
 ERROR_WEIBULL_SHAPE_PARAMETER = "Shape parameter (b) must be positive."
-ERROR_WEIBULL_INITIAL_AMOUNT = (
+ERROR_RELEASABLE_AMOUNT = (
     "Entire releasable amount of drug (M) must be non-negative."
 )
+
+# Hopfenberg model error messages
+ERROR_INVALID_EROSION_CONSTANT = "Erosion rate constant (k0) must be non-negative."
+ERROR_INVALID_INITIAL_RADIUS = "Initial radius or half-thickness (a0) must be positive."
+ERROR_INVALID_GEOMETRY_FACTOR = "Geometry factor (n) must be 1 (slab), 2 (cylinder), or 3 (sphere)."

drux/weibull.py

@@ -4,7 +4,7 @@
 from .base_model import DrugReleaseModel
 from .messages import (
     ERROR_WEIBULL_SCALE_PARAMETER,
-    ERROR_WEIBULL_INITIAL_AMOUNT,
+    ERROR_RELEASABLE_AMOUNT,
     ERROR_WEIBULL_SHAPE_PARAMETER,
 )
 from dataclasses import dataclass
@@ -65,7 +65,7 @@ def _model_function(self, t: float) -> float:
     def _validate_parameters(self) -> None:
         """Validate the parameters of the Weibull model."""
         if self._parameters.M < 0:
-            raise ValueError(ERROR_WEIBULL_INITIAL_AMOUNT)
+            raise ValueError(ERROR_RELEASABLE_AMOUNT)
         if self._parameters.a <= 0:
             raise ValueError(ERROR_WEIBULL_SCALE_PARAMETER)
         if self._parameters.b <= 0:

tests/test_hopfenberg.py

@@ -0,0 +1,151 @@
+"""Tests for the Hopfenberg model implementation in drux package."""
+
+from pytest import raises
+from numpy import isclose
+from re import escape
+from drux import HopfenbergModel
+
+TEST_CASE_NAME = "Hopfenberg model tests"
+M, k0, c0, a0, n = 1, 0.00067, 0.0374, 3.51, 2
+SIM_DURATION, SIM_TIME_STEP = 100, 1
+RELATIVE_TOLERANCE = 1e-1
+
+
+def test_hopfenberg_parameters():
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+    assert model._parameters.M == M
+    assert model._parameters.k0 == k0
+    assert model._parameters.c0 == c0
+    assert model._parameters.a0 == a0
+    assert model._parameters.n == n
+
+
+def test_invalid_parameters():
+    with raises(ValueError, match=escape("Entire releasable amount of drug (M) must be non-negative.")):
+        HopfenbergModel(M=-M, k0=k0, c0=c0, a0=a0, n=n).simulate(duration=SIM_DURATION, time_step=SIM_TIME_STEP)
+
+    with raises(ValueError, match=escape("Erosion rate constant (k0) must be non-negative")):
+        HopfenbergModel(M=M, k0=-k0, c0=c0, a0=a0, n=n).simulate(duration=SIM_DURATION, time_step=SIM_TIME_STEP)
+
+    with raises(ValueError, match=escape("Initial drug concentration (c0) must be positive.")):
+        HopfenbergModel(M=M, k0=k0, c0=-c0, a0=a0, n=n).simulate(duration=SIM_DURATION, time_step=SIM_TIME_STEP)
+
+    with raises(ValueError, match=escape("Initial radius or half-thickness (a0) must be positive.")):
+        HopfenbergModel(M=M, k0=k0, c0=c0, a0=-a0, n=n).simulate(duration=SIM_DURATION, time_step=SIM_TIME_STEP)
+
+    with raises(ValueError, match=escape("Geometry factor (n) must be 1 (slab), 2 (cylinder), or 3 (sphere).")):
+        HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=4).simulate(duration=SIM_DURATION, time_step=SIM_TIME_STEP)
+
+
+def test_repr():
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+    repr_str = repr(model)
+    assert repr_str == f"drux.HopfenbergModel(M={M}, k0={k0}, c0={c0}, a0={a0}, n={n})"
+
+
+def test_hopfenberg_simulation():  # Reference: https://pmc.ncbi.nlm.nih.gov/articles/PMC3500559/
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+    profile = model.simulate(duration=SIM_DURATION, time_step=SIM_TIME_STEP)
+
+    actual_release = [M * (1 - (1 - (k0 * t) / (c0 * a0))**n) for t in range(0, SIM_DURATION+SIM_TIME_STEP, SIM_TIME_STEP)]
+    assert all(isclose(p, r, rtol=RELATIVE_TOLERANCE) for p, r in zip(profile, actual_release))
+
+
+def test_hopfenberg_simulation_errors():
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+
+    with raises(ValueError, match="Duration and time step must be positive values"):
+        model.simulate(duration=-100, time_step=10)
+
+    with raises(ValueError, match="Duration and time step must be positive values"):
+        model.simulate(duration=100, time_step=-10)
+
+    with raises(ValueError, match="Time step cannot be greater than duration"):
+        model.simulate(duration=10, time_step=20)
+
+
+def test_hopfenberg_plot1():
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+    model.simulate(duration=SIM_DURATION, time_step=SIM_TIME_STEP)
+    fig, ax = model.plot()
+    assert fig is not None
+    assert ax is not None
+    assert ax.get_title() == model._plot_parameters["title"]
+    assert ax.get_xlabel() == model._plot_parameters["xlabel"]
+    assert ax.get_ylabel() == model._plot_parameters["ylabel"]
+    assert [text.get_text() for text in ax.get_legend().get_texts()] == [model._plot_parameters["label"]]
+
+
+def test_hopfenberg_plot2():
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+    model.simulate(duration=SIM_DURATION, time_step=SIM_TIME_STEP)
+    fig, ax = model.plot(title="test-title", xlabel="test-xlabel", ylabel="test-ylabel", label="test-label")
+    assert fig is not None
+    assert ax is not None
+    assert ax.get_title() == "test-title"
+    assert ax.get_xlabel() == "test-xlabel"
+    assert ax.get_ylabel() == "test-ylabel"
+    assert [text.get_text() for text in ax.get_legend().get_texts()] == ["test-label"]
+
+
+def test_hopfenberg_plot_error():
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+
+    with raises(ValueError, match=escape("No simulation data available. Run simulate() first.")):
+        model.plot()
+
+    model._time_points = [0]  # manually set time points to simulate error (TODO: it will be caught with prior errors)
+    # manually set a too short profile to simulate error (TODO: it will be caught with prior errors)
+    model._release_profile = [0.0]
+    with raises(ValueError, match="Release profile is too short to calculate release rate."):
+        model.plot()
+
+
+def test_hopfenberg_release_rate():
+    small_timestep = 0.001  # smaller time step for better accuracy in numerical derivative
+    small_duration = 10
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+    model.simulate(duration=small_duration, time_step=small_timestep)
+    rate = model.get_release_rate().tolist()
+    import numpy as np
+    # Ref: https://www.wolframalpha.com/input?i=get+the+derivative+of+%28M+*+%281+-+%281+-+%28k0+*+t%29+%2F+%28c0+*+a0%29%29**n%29%29+with+respect+to+t
+    actual_rate = [
+        k0*M*n*((1-((k0*t)/(a0*c0)))**(n-1))/(a0*c0)
+        for t in np.arange(small_timestep, small_duration + small_timestep, small_timestep)
+    ]  # avoid t=0 to prevent division by zero
+    assert all(isclose(r, ar, rtol=1e-1) for r, ar in zip(rate[2:], actual_rate[1:]))  # skip first two points due to numerical derivative inaccuracies
+
+
+def test_hopfenberg_release_rate_error():
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+
+    with raises(ValueError, match=escape("No simulation data available. Run simulate() first.")):
+        model.get_release_rate()
+
+    model._time_points = [0]  # manually set time points to simulate error (TODO: it will be caught with prior errors)
+    # manually set a too short profile to simulate error (TODO: it will be caught with prior errors)
+    model._release_profile = [0.0]
+    with raises(ValueError, match="Release profile is too short to calculate release rate."):
+        model.get_release_rate()
+
+
+def test_hopfenberg_time_for_release():
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+    model.simulate(duration=SIM_DURATION, time_step=1)
+    tx = model.time_for_release(0.8 * model._release_profile[-1])
+    # Ref: https://www.wolframalpha.com/input?i=solve+for+t+in+%281+-+%281+-+%280.00067*+t%29+%2F+%280.0374+*+3.51%29%29**2%29+%3D+0.8*0.7602750594732341
+    assert isclose(tx, 74, rtol=1e-2)
+
+
+def test_hopfenberg_time_for_release_error():
+    model = HopfenbergModel(M=M, k0=k0, c0=c0, a0=a0, n=n)
+
+    with raises(ValueError, match=escape("No simulation data available. Run simulate() first.")):
+        model.time_for_release(0.5)
+    model.simulate(duration=SIM_DURATION, time_step=SIM_TIME_STEP)
+
+    with raises(ValueError, match="Target release must be non-negative."):
+        model.time_for_release(-0.1)
+
+    with raises(ValueError, match="Target release exceeds maximum release of the simulated duration."):
+        model.time_for_release(model._release_profile[-1] + 0.1)