Условное дозирование в фармакометрике

В этом примере мы покажем, как смоделировать условное дозирование с помощью DiscreteCallbacks. Задача заключается в следующем. В организме пациента действует препарат A(t). Концентрация препарата задана как C(t)=A(t)/V для некоторой константы объема V. В момент времени t=4 пациент обращается в клинику и проходит обследование. Если концентрация препарата в организме ниже 4, пациент получает новую дозу.

Для нашей модели мы будем использовать простое уравнение распада. Мы напишем его в форме на месте, что позволит легко распространить его на более сложные примеры:

using DifferentialEquations
function f(du, u, p, t)
    du[1] = -u[1]
end
u0 = [10.0]
const V = 1
prob = ODEProblem(f, u0, (0.0, 10.0))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 1-element Vector{Float64}:
 10.0

Посмотрим, как выглядит решение без событий.

sol = solve(prob, Tsit5())
using Plots;
gr();
plot(sol)

Мы видим, что в момент времени t=4 пациент должен получить дозу. Закодируем это событие. Необходимо проверить в точке t=4, является ли концентрация u[1]/4 равной <4, и если да, добавить 10 к u[1]. Это можно сделать следующим образом.

condition(u, t, integrator) = t == 4 && u[1] / V < 4
affect!(integrator) = integrator.u[1] += 10
cb = DiscreteCallback(condition, affect!)

DiscreteCallback{typeof(Main.condition), typeof(Main.affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}(Main.condition, Main.affect!, SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 1])

Теперь передадим этот обратный вызов решателю и укажем ему остановиться в момент t=4, чтобы таким образом можно было проверить условие:

sol = solve(prob, Tsit5(), tstops = [4.0], callback = cb)
using Plots;
gr();
plot(sol)

Покажем, что вместо установки значения 10 он действительно добавил 10. Мы могли бы задать значение с помощью affect!(integrator) = integrator.u[1] = 10.

println(sol(4.00000))
println(sol(4.000000000001))

[0.18316389221855156]
[10.183163892208368]

Теперь смоделируем ситуацию с пациентом, у которого скорость распада препарата ниже:

function f(du, u, p, t)
    du[1] = -u[1] / 6
end
u0 = [10.0]
const V = 1
prob = ODEProblem(f, u0, (0.0, 10.0))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 1-element Vector{Float64}:
 10.0

sol = solve(prob, Tsit5())
using Plots;
gr();
plot(sol)

При тех же условиях и при таком же событии этот пациент не получит вторую дозу:

sol = solve(prob, Tsit5(), tstops = [4.0], callback = cb)
using Plots;
gr();
plot(sol)