diff --git a/CHANGELOG.md b/CHANGELOG.md index 574ad21..a680529 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/README.md b/README.md index 87a7aae..61e754b 100644 --- a/README.md +++ b/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) ``` 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) +``` +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
5- S. Dash, "Kinetic modeling on drug release from controlled drug delivery systems," Acta Poloniae Pharmaceutica, 2010.
6- K. H. Ramteke, P. A. Dighe, A. R. Kharat, S. V. Patil, Mathematical models of drug dissolution: A review, Sch. Acad. J. Pharm., vol. 3, no. 5, pp. 388-396, 2014.
7- C. Corsaro, G. Neri, A. M. Mezzasalma, and E. Fazio, "Weibull modeling of controlled drug release from Ag-PMA nanosystems," Polymers, vol. 13, no. 17, p. 2897, 2021.
+
8- H. B. Hopfenberg, "Controlled release from erodible slabs, cylinders, and spheres," in Controlled Release Polymeric Formulations, ACS Symposium Series, vol. 33, pp. 26–32, 1976.
+
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," Research in Pharmaceutical Sciences, vol. 7, no. 1, pp. 57–64, Jan. 2012.
+ ## Show your support ### Star this repo diff --git a/drux/__init__.py b/drux/__init__.py index f2c91ef..e6788d6 100644 --- a/drux/__init__.py +++ b/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 diff --git a/drux/hopfenberg.py b/drux/hopfenberg.py new file mode 100644 index 0000000..89074f9 --- /dev/null +++ b/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) diff --git a/drux/messages.py b/drux/messages.py index 1ca0b1d..9d7b189 100644 --- a/drux/messages.py +++ b/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)." diff --git a/drux/weibull.py b/drux/weibull.py index 522ec8a..8754df2 100644 --- a/drux/weibull.py +++ b/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: diff --git a/otherfiles/hopfenberg_plot.png b/otherfiles/hopfenberg_plot.png new file mode 100644 index 0000000..dd9825f Binary files /dev/null and b/otherfiles/hopfenberg_plot.png differ diff --git a/tests/test_hopfenberg.py b/tests/test_hopfenberg.py new file mode 100644 index 0000000..488fbc5 --- /dev/null +++ b/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)