|
| 1 | +# |
| 2 | +# Stochastic logistic model. |
| 3 | +# |
| 4 | +# This file is part of PINTS (https://github.com/pints-team/pints/) which is |
| 5 | +# released under the BSD 3-clause license. See accompanying LICENSE.md for |
| 6 | +# copyright notice and full license details. |
| 7 | +# |
| 8 | +from __future__ import absolute_import, division |
| 9 | +from __future__ import print_function, unicode_literals |
| 10 | +import numpy as np |
| 11 | +from scipy.interpolate import interp1d |
| 12 | +import pints |
| 13 | + |
| 14 | +from . import ToyModel |
| 15 | + |
| 16 | + |
| 17 | +class StochasticLogisticModel(pints.ForwardModel, ToyModel): |
| 18 | + r""" |
| 19 | + This model describes the growth of a population of individuals, where the |
| 20 | + birth rate per capita, initially :math:`b_0`, decreases to :math:`0` as the |
| 21 | + population size, :math:`\mathcal{C}(t)`, starting from an initial |
| 22 | + population size, :math:`n_0`, approaches a carrying capacity, :math:`k`. |
| 23 | + This process follows a rate according to [1]_ |
| 24 | +
|
| 25 | + .. math:: |
| 26 | + A \xrightarrow{b_0(1-\frac{\mathcal{C}(t)}{k})} 2A. |
| 27 | +
|
| 28 | + The model is simulated using the Gillespie stochastic simulation algorithm |
| 29 | + [2]_, [3]_. |
| 30 | +
|
| 31 | + *Extends:* :class:`pints.ForwardModel`, :class:`pints.toy.ToyModel`. |
| 32 | +
|
| 33 | + Parameters |
| 34 | + ---------- |
| 35 | + initial_molecule_count : float |
| 36 | + Sets the initial population size :math:`n_0`. |
| 37 | +
|
| 38 | + References |
| 39 | + ---------- |
| 40 | + .. [1] Simpson, M. et al. 2019. Process noise distinguishes between |
| 41 | + indistinguishable population dynamics. bioRxiv. |
| 42 | + https://doi.org/10.1101/533182 |
| 43 | + .. [2] Gillespie, D. 1976. A General Method for Numerically Simulating the |
| 44 | + Stochastic Time Evolution of Coupled Chemical Reactions. |
| 45 | + Journal of Computational Physics. 22 (4): 403-434. |
| 46 | + https://doi.org/10.1016/0021-9991(76)90041-3 |
| 47 | + .. [3] Erban R. et al. 2007. A practical guide to stochastic simulations |
| 48 | + of reaction-diffusion processes. arXiv. |
| 49 | + https://arxiv.org/abs/0704.1908v2 |
| 50 | + """ |
| 51 | + |
| 52 | + def __init__(self, initial_molecule_count=50): |
| 53 | + super(StochasticLogisticModel, self).__init__() |
| 54 | + self._n0 = float(initial_molecule_count) |
| 55 | + if self._n0 < 0: |
| 56 | + raise ValueError('Initial molecule count cannot be negative.') |
| 57 | + |
| 58 | + def n_parameters(self): |
| 59 | + """ See :meth:`pints.ForwardModel.n_parameters()`. """ |
| 60 | + return 2 |
| 61 | + |
| 62 | + def _simulate_raw(self, parameters): |
| 63 | + """ |
| 64 | + Returns tuple (raw times, population sizes) when reactions occur. |
| 65 | + """ |
| 66 | + parameters = np.asarray(parameters) |
| 67 | + if len(parameters) != self.n_parameters(): |
| 68 | + raise ValueError('This model should have only 2 parameters.') |
| 69 | + b = parameters[0] |
| 70 | + k = parameters[1] |
| 71 | + if b <= 0: |
| 72 | + raise ValueError('Rate constant must be positive.') |
| 73 | + |
| 74 | + # Initial time and count |
| 75 | + t = 0 |
| 76 | + a = self._n0 |
| 77 | + |
| 78 | + # Run stochastic logistic birth-only algorithm, calculating time until |
| 79 | + # next reaction and increasing population count by 1 at that time |
| 80 | + mol_count = [a] |
| 81 | + time = [t] |
| 82 | + while a < k: |
| 83 | + r = np.random.uniform(0, 1) |
| 84 | + t += np.log(1 / r) / (a * b * (1 - a / k)) |
| 85 | + a = a + 1 |
| 86 | + time.append(t) |
| 87 | + mol_count.append(a) |
| 88 | + return time, mol_count |
| 89 | + |
| 90 | + def _interpolate_values(self, time, pop_size, output_times, parameters): |
| 91 | + """ |
| 92 | + Takes raw times and population size values as inputs and outputs |
| 93 | + interpolated values at output_times. |
| 94 | + """ |
| 95 | + # Interpolate as step function, increasing pop_size by 1 at each |
| 96 | + # event time point |
| 97 | + interp_func = interp1d(time, pop_size, kind='previous') |
| 98 | + |
| 99 | + # Compute population size values at given time points using f1 |
| 100 | + # at any time beyond the last event, pop_size = k |
| 101 | + values = interp_func(output_times[np.where(output_times <= time[-1])]) |
| 102 | + zero_vector = np.full( |
| 103 | + len(output_times[np.where(output_times > time[-1])]), |
| 104 | + parameters[1]) |
| 105 | + values = np.concatenate((values, zero_vector)) |
| 106 | + return values |
| 107 | + |
| 108 | + def simulate(self, parameters, times): |
| 109 | + """ See :meth:`pints.ForwardModel.simulate()`. """ |
| 110 | + times = np.asarray(times) |
| 111 | + if np.any(times < 0): |
| 112 | + raise ValueError('Negative times are not allowed.') |
| 113 | + if self._n0 == 0: |
| 114 | + return np.zeros(times.shape) |
| 115 | + |
| 116 | + # run Gillespie |
| 117 | + time, pop_size = self._simulate_raw(parameters) |
| 118 | + |
| 119 | + # interpolate |
| 120 | + values = self._interpolate_values(time, pop_size, times, parameters) |
| 121 | + return values |
| 122 | + |
| 123 | + def mean(self, parameters, times): |
| 124 | + r""" |
| 125 | + Computes the deterministic mean of infinitely many stochastic |
| 126 | + simulations with times :math:`t` and parameters (:math:`b`, :math:`k`), |
| 127 | + which follows: |
| 128 | + :math:`\frac{kC(0)}{C(0) + (k - C(0)) \exp(-bt)}`. |
| 129 | +
|
| 130 | + Returns an array with the same length as `times`. |
| 131 | + """ |
| 132 | + parameters = np.asarray(parameters) |
| 133 | + if len(parameters) != self.n_parameters(): |
| 134 | + raise ValueError('This model should have only 2 parameters.') |
| 135 | + |
| 136 | + b = parameters[0] |
| 137 | + if b <= 0: |
| 138 | + raise ValueError('Rate constant must be positive.') |
| 139 | + |
| 140 | + k = parameters[1] |
| 141 | + if k <= 0: |
| 142 | + raise ValueError("Carrying capacity must be positive") |
| 143 | + |
| 144 | + times = np.asarray(times) |
| 145 | + if np.any(times < 0): |
| 146 | + raise ValueError('Negative times are not allowed.') |
| 147 | + c0 = self._n0 |
| 148 | + return (c0 * k) / (c0 + np.exp(-b * times) * (k - c0)) |
| 149 | + |
| 150 | + def variance(self, parameters, times): |
| 151 | + r""" |
| 152 | + Returns the deterministic variance of infinitely many stochastic |
| 153 | + simulations. |
| 154 | + """ |
| 155 | + raise NotImplementedError |
| 156 | + |
| 157 | + def suggested_parameters(self): |
| 158 | + """ See :meth:`pints.toy.ToyModel.suggested_parameters()`. """ |
| 159 | + return np.array([0.1, 500]) |
| 160 | + |
| 161 | + def suggested_times(self): |
| 162 | + """ See :meth:`pints.toy.ToyModel.suggested_times()`.""" |
| 163 | + return np.linspace(0, 100, 101) |
0 commit comments