How to force my differential equations give a bounded solution?

Question

I have modeled the interaction of two physical quantities, $S$ and $B$, by the following differential equations (the second one is a delay differential equation):

$$S'(t) = 0.31 S(t) \Big( 1 - \frac{S(t)}{6.19} \Big)- \frac{a}{1 + B(t)} S(t),$$

$$B'(t) = c (1 - B(t)) \theta(18 - t) - d B(t) S(t) - 3 c (1 - B(t - 18)),$$

where $a$, $c$, and $d$ are some constants to be determined from the experimental data and $\theta$ is the Heaviside function. The parameters are found to be: $a = 0.72$, $c = 0.19$, and $d = 0.24$; and also, the initial conditions read: $S(0) = 1.25$ and $B(0) = 1.14$.

The solutions (plotted with Mathematica) are:

Unprotect[HeavisideTheta]; HeavisideTheta[0.] = 0.5; Protect[HeavisideTheta];
so = NDSolve[{ss'[t] + 0.72/(1 + bb[t]) ss[t] - 0.31 ss[t] (1 - ss[t]/6.19) == 0, 
bb'[t] - 0.19 (1 - bb[t]) HeavisideTheta[18 - t] + 0.24 bb[t] ss[t] + 0.57 (1 - bb[t - 18]) == 0, 
ss[t /; t < 0] == 1.25, bb[t /; t < 0] == 1.14}, {ss, bb}, {t, 0, 24}];

Plot[Evaluate[bb[t] /. so], {t, 0, 24}, PlotRange -> All, AxesOrigin -> {0, 0}, PlotLabel -> "B(t)"]
Plot[Evaluate[ss[t] /. so], {t, 0, 24}, PlotRange -> All, AxesOrigin -> {0, 0}, PlotLabel -> "S(t)"]

Image 1

Everything looks good, since I'm only interested in the time interval $0 < t < 24$ and the curves correctly predict the behavior of my system (experiment). However, if we simulate the differential equations further and for a larger time interval, let's say, till $t = 30$, the solutions look like:

Image 2

This is not desirable as $B$ is a physical quantity and cannot be negative.

My question is:

How can I introduce a cut-off in my differential equations to force the $B$ solution do not cross the axis and to be bounded and remain zero at the end?

Any help is appreciated.

Answer 1

You could just modify your second equation once $t>24$, e.g., $$ B'(t) = \theta(24-t) [c (1 - B(t)) \theta(18 - t) - d B(t) S(t) - 3 c (1 - B(t - 18))]-\frac{1}{2}\theta(t-24)B(t) $$ If a smoother modification is desired, one could, e.g., replace the newly introduced step functions as $\theta(x) \rightarrow (\tanh(10x)+1)/2$.