<jupyter_start><jupyter_text>Identifiability of Structural Growth Model in Absence of CompoundThe tumour growth inhibiting effects of Erlotinib and Gefitinib were modelled with a population PKPD model in [1]. A population PKPD model is a hierarchical model which consists of a structural model, a population model, and an error model. Each sub-model captures a different aspect of the tumour growth inhibition biology, and is parametrised by a set of parameters. In this notebook we explore the identifiability of the structural model reported in [1] in absence of the compound. In particular we investigate the effect of log-transforms and non-dimensionalisation of the model parameters on the stability of the inference.   Structural growth model in absence of the compoundIn a PKPD model the structural model is a mechanistic or empirical description of the drug's pharmacokinetics (PK) and pharmacodynamics (PD). In particular, the PD model needs to capture not only the effects of the drug on a specfic biomarker or pathology read-out, but also model the biology in absence of the compund. In this study the pathology read-out is the tumour size. In [1] a structural model is proposed, which assumes that in absence of the compound the tumour grows exponentially for small tumour volumes, and linearly for large tumour volumes\begin{equation*}    \frac{\text{d}V^s_T}{\text{d}t} = \frac{2\lambda _0\lambda _1 V^s_T}{2\lambda _0 V^s_T + \lambda _1}.\end{equation*}Here, $V^s_T$ is the predicted tumour volume by the structural model, $\lambda _0$ is the exponential growth rate, and $\lambda _1$ is the linear growth rate. The tumour growth of an individual in absence of the drug is thus parameterised by three parameters\begin{equation*}    \psi = (V_0, \lambda _0, \lambda _1),\end{equation*}where $V_0$ is the initial tumour volume, $V^s_T(t=0, \psi) = V_0$.  IdentifiabilityThe identifiability of parametric models is a multifacetted problem and has been addressed many times in the literature, e.g. [2, 3, 4]. Generally, identifiability of a model refers to the ability to infer its parameters from data. In particular, under the assumption that the model is able to generate the data with a given set of model parameters, identifiability addresses the question of whether those parameters can be recovered from the data. In the PKPD modelling context, the identifiability of the structural model is of particular importance, as model parameters are often biologically interpretable. This allows to characterise biochemical properties of the compound based on the inferred model parameters, and even to translate PKPD models from preclinical to clinical application. In this study, we are interested to learn the structural model parameters $\psi $ of the tumour growth in absence of the compound from two *in vivo* experiments reported in [1]. Here, patient-derived tumour explants LXF A677 (adenocarcinoma of the lung) and cell line-derived tumour xenografts VXF A431 (vulva cancer) were implanted in mice. The tumour growth was monitored over a period of 30 days, see the [Control Growth Data notebook](https://nbviewer.jupyter.org/github/DavAug/ErlotinibGefitinib/blob/master/notebooks/control_growth/data_preparation.ipynb) for details.<jupyter_code>#
# Visualise control growth data.
#

import os

import pandas as pd
import plotly.graph_objects as go


# Create figure
fig = go.Figure()

# Import data
# Get path of current working directory
path = os.getcwd()

# Import LXF A677 control growth data
lxf_data = pd.read_csv(path + '/data/lxf_control_growth.csv')

# Import VXF A341 control growth data
vxf_data = pd.read_csv(path + '/data/vxf_control_growth.csv')

# Scatter plot LXF A677 time-series data for each mouse
mouse_ids = lxf_data['#ID'].unique()
for id_m in mouse_ids:
    # Create mask for mouse
    mask = lxf_data['#ID'] == id_m

    # Get time points for mouse
    times = lxf_data['TIME in day'][mask]

    # Get observed tumour volumes for mouse
    observed_volumes = lxf_data['TUMOUR VOLUME in cm^3'][mask]

    # Plot data
    fig.add_trace(go.Scatter(
        x=times,
        y=observed_volumes,
        legendgroup="LXF A677",
        name="LXF A677",
        showlegend=True if id_m == mouse_ids[0] else False,
        hovertemplate=
            "<b>%s ID: %d</b><br>" % ("LXF A677", id_m) +
            "Time: %{x:} day<br>" +
            "Tumour volume: %{y:.02f} cm^3<br>" +
            "<extra></extra>",
        mode="markers",
        marker=dict(
            symbol='circle',
            opacity=0.7,
            line=dict(color='black', width=1))
    ))

# Scatter plot VXF A341 time-series data for each mouse
mouse_ids = vxf_data['#ID'].unique()
for id_m in mouse_ids:
    # Create mask for mouse
    mask = vxf_data['#ID'] == id_m

    # Get time points for mouse
    times = vxf_data['TIME in day'][mask]

    # Get observed tumour volumes for mouse
    observed_volumes = vxf_data['TUMOUR VOLUME in cm^3'][mask]

    # Plot data
    fig.add_trace(go.Scatter(
        x=times,
        y=observed_volumes,
        legendgroup="VXF A341",
        name="VXF A341",
        showlegend=True if id_m == mouse_ids[0] else False,
        hovertemplate=
            "<b>%s ID: %d</b><br>" % ("VXF A341", id_m) +
            "Time: %{x:} day<br>" +
            "Tumour volume: %{y:.02f} cm^3<br>" +
            "<extra></extra>",
        mode="markers",
        marker=dict(
            symbol='star',
            opacity=0.7, 
            line=dict(color='black', width=1))
    ))

# Set X, Y axis and figure size
fig.update_layout(
    autosize=True,
    xaxis_title=r'$\text{Time in day}$',
    yaxis_title=r'$\text{Tumour volume in cm}^3$',
    template="plotly_white")

# Add switch between linear and log y-scale
fig.update_layout(
    updatemenus=[
        dict(
            type = "buttons",
            direction = "left",
            buttons=list([
                dict(
                    args=[{"yaxis.type": "linear"}],
                    label="Linear y-scale",
                    method="relayout"
                ),
                dict(
                    args=[{"yaxis.type": "log"}],
                    label="Log y-scale",
                    method="relayout"
                )
            ]),
            pad={"r": 0, "t": -10},
            showactive=True,
            x=0.0,
            xanchor="left",
            y=1.15,
            yanchor="top"
        ),
    ]
)

# Show figure
fig.show()<jupyter_output><empty_output><jupyter_text>**Figure 1:** Untreated tumour growth of patient-derived tumour explants LXF A677 (adenocarcinoma of the lung) and cell line-derived tumour xenografts VXF A431 (vulva cancer) in mice. The colouring of the data points indicates that the measurements belong to the same mouse. Mouse ID and further information can be explored by hovering over the data points. Naive check: fit strucurla model with SD to each indivual separartelt<jupyter_code>#
# Create structural model.
#

import pkpd.model as m


# Create model
model = m.create_tumour_growth_model()

# Show model
print(model.code())
#
# Define pints model wrapper such that myokit model can be used for inference.
#

import myokit
import pints

from pkpd import model as model


# Wrap myokit model, so it can be used with pints
class PintsModel(pints.ForwardModel):
    def __init__(self):
        # Create myokit model
        model = m.create_tumour_growth_model()

        # Create simulator
        self.sim = myokit.Simulation(model)

    def n_parameters(self):
        """
        Number of parameters to fit. Here initial V^s_T, lambda_0, lambda_1
        """
        return 3

    def n_outputs(self):
        return 1

    def simulate(self, parameters, times):
        # Reset simulation
        self.sim.reset()

        # Sort input parameters
        initial_volume, lambda_0, lambda_1 = parameters

        # Set initial condition
        self.sim.set_state([initial_volume])

        # Set growth constants
        self.sim.set_constant('central.lambda_0', lambda_0)
        self.sim.set_constant('central.lambda_1', lambda_1)

        # Define logged variable
        loggedVariable = 'central.volume_t'

        # Simulate
        output = self.sim.run(times[-1] + 1, log=[loggedVariable], log_times=times)
        result = output[loggedVariable]

        return np.array(result)
#
# Naive attempt to find an optimal parameter set psi for each mouse in the LXF A677 and VXF A431 population.
#
# This cell needs above defined wrapped myokit model:
# [PintsModel]
#

import os

import numpy as np
import pandas as pd
import pints


# Import data
# Get path of current working directory
path = os.getcwd()

# Export cleaned LXF A677 control growth data
lxf_data = pd.read_csv(path + '/data/lxf_control_growth.csv')
n_mice_lxf = len(lxf_data['#ID'].unique())

# Export cleaned VXF A341 control growth data
vxf_data = pd.read_csv(path + '/data/vxf_control_growth.csv')
n_mice_vxf = len(vxf_data['#ID'].unique())

# Define container for the structural model estimates
# Shape (n_mice, n_parameters)
n_parameters = 3
lxf_mouse_parameters = np.empty(shape=(n_mice_lxf, n_parameters))
vxf_mouse_parameters = np.empty(shape=(n_mice_vxf, n_parameters)) 

# Define arbitrary starting point for the optimisations
initial_parameters = [1, 1, 1]

# Find mouse parameters for LXF A677 population
mouse_ids = lxf_data['#ID'].unique()
for index, mouse_id in enumerate(mouse_ids):
    # Create mask for mouse with specfied ID
    mouse_mask = lxf_data['#ID'] == mouse_id

    # Get relevant time points
    times = lxf_data[mouse_mask]['TIME in day'].to_numpy()

    # Get measured tumour volumes
    observed_volumes = lxf_data[mouse_mask]['TUMOUR VOLUME in cm^3'].to_numpy()

    # Create inverse problem
    problem = pints.SingleOutputProblem(PintsModel(), times, observed_volumes)

    # Create sum of squares error objective function
    error = pints.SumOfSquaresError(problem)

    # Create optimisation controller with a CMA-ES optimiser
    optimiser = pints.OptimisationController(
        function=error,
        x0=initial_parameters,
        method=pints.CMAES)

    # Disable logging mode
    optimiser.set_log_to_screen(False)

    # Parallelise optimisation
    optimiser.set_parallel(True)

    # Find optimal parameters
    estimates, _ = optimiser.run()

    # Save estimates
    lxf_mouse_parameters[index, :] = estimates
#
# Solve structural model for inferred model parameters.
#
# This cell needs the above defined wrapped myokit model, and the inferred model parameters:
# [PintsModel, lxf_mouse_parameters, vxf_mouse_parameters]
#

import numpy as np


# Create tumour growth model
model = PintsModel()

# Define simulation time points in day
times = np.linspace(start=0, stop=30, num=200)
n_times = len(times)

# Create container for simulated tumour growth
# Shape (n_mice, n_times)
n_mice = len(lxf_mouse_parameters)
lxf_tumour_growth = np.empty(shape=(n_mice, n_times))

# Solve structural model for LXF A677 population
for mouse_id, mouse_params in enumerate(lxf_mouse_parameters):
    # Simulate mouse tumour growth
    lxf_tumour_growth[mouse_id, :] = model.simulate(parameters=mouse_params, times=times)

#
# Visualise simulated tumour growth together with control growth data.
#
# This cell needs the above simulated tumour growth and the control growth data:
# [times, lxf_tumour_growth, lxf_data]
#

import pandas as pd
import plotly.colors
import plotly.graph_objects as go


# Get number of individual mice
n_mice_lxf = len(lxf_tumour_growth)

# Define colorscheme
colors = plotly.colors.qualitative.Plotly[:n_mice_lxf]

# Create figure
fig = go.Figure()

# Scatter plot LXF A677 time-series data for each mouse
mouse_ids = lxf_data['#ID'].unique()
for index, id_m in enumerate(mouse_ids):
    # Create mask for mouse
    mask = lxf_data['#ID'] == id_m

    # Get time points for mouse
    observed_times = lxf_data['TIME in day'][mask]

    # Get observed tumour volumes for mouse
    observed_volumes = lxf_data['TUMOUR VOLUME in cm^3'][mask]

    # Get simulated tumour volumes for mouse
    simulated_volumes = lxf_tumour_growth[index, :]

    # Plot data
    fig.add_trace(go.Scatter(
        x=observed_times,
        y=observed_volumes,
        legendgroup="ID: %d" % id_m,
        name="ID: %d" % id_m,
        showlegend=True,
        hovertemplate=
            "<b>Measurements %s</b><br>" % ("LXF A677") +
            "ID: %d<br>" % id_m +
            "Time: %{x:} day<br>" +
            "Tumour volume: %{y:.02f} cm^3<br>" +
            "<extra></extra>",
        mode="markers",
        marker=dict(
            symbol='circle',
            opacity=0.7,
            line=dict(color='black', width=1),
            color=colors[index])
    ))

    # Plot simulated growth
    fig.add_trace(go.Scatter(
        x=times,
        y=simulated_volumes,
        legendgroup="ID: %d" % id_m,
        name="ID: %d" % id_m,
        showlegend=False,
        hovertemplate=
            "<b>Simulation %s</b><br>" % ("LXF A677") +
            "ID: %d<br>" % id_m +
            "Time: %{x:.0f} day<br>" +
            "Tumour volume: %{y:.02f} cm^3<br>" +
            "<extra></extra>",
        mode="lines",
        line=dict(color=colors[index])
    ))

# # Scatter plot VXF A341 time-series data for each mouse
# mouse_ids = vxf_data['#ID'].unique()
# for id_m in mouse_ids:
#     # Create mask for mouse
#     mask = vxf_data['#ID'] == id_m

#     # Get time points for mouse
#     times = vxf_data['TIME in day'][mask]

#     # Get observed tumour volumes for mouse
#     observed_volumes = vxf_data['TUMOUR VOLUME in cm^3'][mask]

#     # Plot data
#     fig.add_trace(go.Scatter(
#         x=times,
#         y=observed_volumes,
#         legendgroup="VXF A341",
#         name="VXF A341",
#         showlegend=True if id_m == mouse_ids[0] else False,
#         hovertemplate=
#             "<b>%s ID: %d</b><br>" % ("VXF A341", id_m) +
#             "Time: %{x:} day<br>" +
#             "Tumour volume: %{y:.02f} cm^3<br>" +
#             "<extra></extra>",
#         mode="markers",
#         marker=dict(
#             symbol='star',
#             opacity=0.7, 
#             line=dict(color='black', width=1))
#     ))

# Set X, Y axis and figure size
fig.update_layout(
    autosize=True,
    xaxis_title=r'$\text{Time in day}$',
    yaxis_title=r'$\text{Tumour volume in cm}^3$',
    template="plotly_white")

# Add switch between linear and log y-scale
fig.update_layout(
    updatemenus=[
        dict(
            type = "buttons",
            direction = "left",
            buttons=list([
                dict(
                    args=[{"yaxis.type": "linear"}],
                    label="Linear y-scale",
                    method="relayout"
                ),
                dict(
                    args=[{"yaxis.type": "log"}],
                    label="Log y-scale",
                    method="relayout"
                )
            ]),
            pad={"r": 0, "t": -10},
            showactive=True,
            x=0.0,
            xanchor="left",
            y=1.15,
            yanchor="top"
        ),
    ]
)

# Show figure
fig.show()<jupyter_output><empty_output><jupyter_text>Causes for non-identifiabilityIn theory there are only two scenarios in which a structural model may be rendered non-identifiable. On the one hand, the model may be overparametrised, in the sense that two distinct sets of model parameters $\psi $ and $\psi '$ generate an identical model output, $V^s_T(t, \psi) = V^s_T(t, \psi ')$, for $\psi \neq \psi '$. If that is the case, a hypothetical data-generating set of model parameters cannot distingished from another set of model parameters, which happens to produce the same model output. In that case the model is said to be structurally non-identifiable [2]. On the other hand, the data may simply not be rich enough to distingish two distinct parameter sets $\psi $ and $\psi '$. In other words, the model output for all time and sampled time points.Numerical issues. There are two ways in which an ODE model may be rendered non-identifiable: a model can be structurally non-identifiable; and it can be practically non-identifiable [5].We are interested in the parameters $V_0$, $\lambda _0$ and $\lambda _1$ that best describe the observed tumour growth in untreated mice, see the dedicated [notebook](https://nbviewer.jupyter.org/github/DavAug/ErlotinibGefitinib/blob/master/notebooks/control_growth/data_preparation.ipynb) for details on the data. Here $V_0$ is the initial tumour volume at time $t=0$ according to the structural model, i.e. $V^s_T(t=0)=V_0$. In order to find good parameters we will use Bayesian inference. Bayesian inference constructs distributions for the parameters $V_0$, $\lambda _0$ and $\lambda _1$, indicating the probability with which certain parameter values have generated the data.In this notebook we want to make sure that - identifiability: technical challenges versus theoretical problems- discuss reasons for non-dimensionalisation    - in theory Bayesian inference works for any model    - in practice performance is a lot improved by by having parameters of equal scale (most optimisation have only one step size)    - avoid hard boundaries, discuss why- model can be such that it's overparametrised- or data does not provide enough information about data  Implementation of naïve pooled modelWe are using [myokit](http://myokit.org/) for the implementation of the structural model. Myokit enables us to solve the structural model ODE with an adaptive numerical solver called CVODE [3]. To implement the error model and perform the inference we are using [pints](https://pints.readthedocs.io/). Note that in general the quality of the inference of $\psi $ and $\theta _V$ can be significantly improved when all parameters are appropriately transformed. We will however choose to not transform the parameetrs at first, to illustrate how the inference may be stabilised with transformations.  Implementation of structural model in myokitWe have implemented the structural model in myokit with untransformed parameters in a separate [module](https://github.com/DavAug/ErlotinibGefitinib/blob/master/pkpd/model.py). The structural model can now be created by calling ```pkpd.model.create_tumour_growth_model()```.<jupyter_code>#
# Implementing the structural model in myokit.
#

from pkpd import model as m


# Create model
model = m.create_tumour_growth_model()

# Print structural model
print(model.code())<jupyter_output>[[model]]
# Initial values
central.volume_t = 0

[central]
lambda_0 = 0
    in [1/day]
lambda_1 = 1
    in [cm^3/day]
time = 0 bind time
    in [day]
dot(volume_t) = 2 * (lambda_0 * (lambda_1 * volume_t)) / (2 * (lambda_0 * volume_t) + lambda_1)
    in [cm^3]<jupyter_text>Implementation of error model We use pints to implement the error model indirectly by defining a log-likelihood for the model parameters $\psi $ and $\theta _V$\begin{equation*}L(\psi , \theta _V | V^{\text{obs}}_T) = \mathbb{P}(V_T | \psi , \theta _V) \Big | _{V_T = V^{\text{obs}}_T} .\end{equation*}Note that here the likelihood is mathematically equivalent to the above defined distribution of tumour growth curves.To use the structural model in pints we first need to wrap the myokit model by a `pints.ForwardModel`.<jupyter_code>#
# Define pints model wrapper such that myokit model can be used for inference.
#

import myokit
import pints

from pkpd import model as model


# Wrap myokit model, so it can be used with pints
class PintsModel(pints.ForwardModel):
    def __init__(self):
        # Create myokit model
        model = m.create_tumour_growth_model()

        # Create simulator
        self.sim = myokit.Simulation(model)

    def n_parameters(self):
        """
        Number of parameters to fit. Here initial V^s_T, lambda_0, lambda_1
        """
        return 3

    def n_outputs(self):
        return 1

    def simulate(self, parameters, times):
        # Reset simulation
        self.sim.reset()

        # Sort input parameters
        initial_volume, lambda_0, lambda_1 = parameters

        # Set initial condition
        self.sim.set_state([initial_volume])

        # Set growth constants
        self.sim.set_constant('central.lambda_0', lambda_0)
        self.sim.set_constant('central.lambda_1', lambda_1)

        # Define logged variable
        loggedVariable = 'central.volume_t'

        # Simulate
        output = self.sim.run(times[-1] + 1, log=[loggedVariable], log_times=times)
        result = output[loggedVariable]

        return np.array(result)<jupyter_output><empty_output><jupyter_text>We also need to import the tumour growth data that we cleaned in a previous [notebook](https://github.com/DavAug/ErlotinibGefitinib/blob/master/notebooks/control_growth/data_preparation.ipynb) to define the likelihood of the model parameters.<jupyter_code>#
# Load cleaned LXF and VXF data sets.
#

import pandas as pd


# Get path of current working directory
path = os.getcwd()

# Read cleaned LXF A677 control growth data
lxf_data = pd.read_csv(path + '/data/lxf_control_growth.csv')

# Read cleaned VXF A341 control growth data
vxf_data = pd.read_csv(path + '/data/vxf_control_growth.csv')<jupyter_output><empty_output><jupyter_text>With a `pints.ForwardModel` implementing the structural model for the tumour growth and the time-series data of the control tumour growth, we can now construct the likelihood for the model parameters $\psi $ and $\theta $.(Note that we still have to create likelihoods for the individual mice, because myokit does not allow multiple measurements at the same time point.)  Likelihood for LXF A677 tumour growth model parameters<jupyter_code>#
# Construct likelihood for LXF A677 tumour growth model parameters.
#


import numpy as np
import pandas as pd
import pints


# Create inverse problem for each mouse ID
mouse_ids = lxf_data['#ID'].unique()
log_likelihoods = []
for ids in mouse_ids:
    # Create mask for mouse with specfied ID
    mouse_mask = lxf_data['#ID'] == ids

    # Get relevant time points
    times = lxf_data[mouse_mask]['TIME in day'].to_numpy()

    # Get measured tumour volumes
    observed_volumes = lxf_data[mouse_mask]['TUMOUR VOLUME in cm^3'].to_numpy()

    # Create inverse problem
    problem = pints.SingleOutputProblem(PintsModel(), times, observed_volumes)

    # Create Gaussian log-likelihood TODO: Change to combined error model
    log_likelihoods.append(pints.GaussianLogLikelihood(problem))

# Create one log_likelihood for the inference from the individual problems
lxf_log_likelihood = pints.SumOfIndependentLogPDFs(log_likelihoods)<jupyter_output><empty_output><jupyter_text>Likelihood for VXF A341 tumour growth model parameters<jupyter_code>#
# Construct likelihood for VXF A341 tumour growth model parameters.
#


import numpy as np
import pandas as pd
import pints


# Create inverse problem for each mouse ID
mouse_ids = vxf_data['#ID'].unique()
log_likelihoods = []
for ids in mouse_ids:
    # Create mask for mouse with specfied ID
    mouse_mask = vxf_data['#ID'] == ids

    # Get relevant time points
    times = vxf_data[mouse_mask]['TIME in day'].to_numpy()

    # Get measured tumour volumes
    observed_volumes = vxf_data[mouse_mask]['TUMOUR VOLUME in cm^3'].to_numpy()

    # Create inverse problem
    problem = pints.SingleOutputProblem(PintsModel(), times, observed_volumes)

    # Create Gaussian log-likelihood TODO: Change to combined error model
    log_likelihoods.append(pints.GaussianLogLikelihood(problem))

# Create one log_likelihood for the inference from the individual problems
vxf_log_likelihood = pints.SumOfIndependentLogPDFs(log_likelihoods)<jupyter_output><empty_output><jupyter_text>Identifiability checksome explanation<jupyter_code>#
# Check identifiability of LXF A677 problem.
#

import pints


# Get number of optimisation runs and number of parameters
n = 5
n_parameters = lxf_log_likelihood.n_parameters()

# Initial guess of parameters
parameters = np.array([
    [10E-2, 10E-2, 10E-2, 10E-2],
    [10E-1, 10E-1, 10E-1, 10E-1],
    [10E0, 10E0, 10E0, 10E0],
    [10E1, 10E1, 10E1, 10E1],
    [10E2, 10E2, 10E2, 10E2]])

# Create estimate container
estimates = np.empty(shape=(n, n_parameters))

print('Running...')
# estimate parameters
for i in range(n):
    print('Round %d' % (i+1))

    # Initial guess of parameters
    params = parameters[i, :]

    # Choose optimisation method
    optimiser = pints.CMAES

    # Create optimisation object
    opt = pints.OptimisationController(
        function=lxf_log_likelihood,
        x0=params,
        method=optimiser)

    # Disable logging mode
    opt.set_log_to_screen(False)

    # Parallelise optimisation
    opt.set_parallel(True)

    # Optimise likelihood
    est, _ = opt.run()

    # Save estimates
    estimates[i, :] = est

print('Done!')

print(' ')
print('Estimates: ')
print('Initial tumour volume [cm^3]: ', estimates[:, 0])
print('Exponential growth rate \lambda _0 [1/day]: ', estimates[:, 1])
print('Linear growth rate \lambda _1 [cm^3/day]: ', estimates[:, 2])
print('Standard deviation of base-level noise [cm^3]: ', estimates[:, 3])<jupyter_output>Running...
Round 1
Round 2
Round 3
Round 4
Round 5
Done!
 
Estimates: 
Initial tumour volume [cm^3]:  [1.34796622 0.93024172 1.34753463 1.34111851 1.27354206]
Exponential growth rate \lambda _0 [1/day]:  [ 3.66072941e+02 -2.50009666e-01  1.82867188e+02 -1.01997319e+01
 -3.28346683e-02]
Linear growth rate \lambda _1 [cm^3/day]:  [  0.3587859    0.31194863   0.35883914   0.35750309 144.84249923]
Standard-deviation of base-level noise [cm^3]:  [4.51403325 4.50053914 4.51405733 4.51374813 8.34865767]<jupyter_text>Log-transform structural model parameters<jupyter_code>#
# Define pints model wrapper with log-transformed parameters.
#

import myokit
import numpy as np
import pints

from pkpd import model as model


# Wrap myokit model, so it can be used with pints
class PintsModel(pints.ForwardModel):
    def __init__(self):
        # Create myokit model
        model = m.create_tumour_growth_model()

        # Create simulator
        self.sim = myokit.Simulation(model)

    def n_parameters(self):
        """
        Number of parameters to fit. Here initial V^s_T, lambda_0, lambda_1
        """
        return 3

    def n_outputs(self):
        return 1

    def simulate(self, log_parameters, times):
        # Reset simulation
        self.sim.reset()

        # Sort input parameters
        initial_volume, lambda_0, lambda_1 = np.exp(log_parameters)

        # Set initial condition
        self.sim.set_state([initial_volume])

        # Set growth constants
        self.sim.set_constant('central.lambda_0', lambda_0)
        self.sim.set_constant('central.lambda_1', lambda_1)

        # Define logged variable
        loggedVariable = 'central.volume_t'

        # Simulate
        output = self.sim.run(times[-1] + 1, log=[loggedVariable], log_times=times)
        result = output[loggedVariable]

        return np.array(result)
#
# Construct likelihood for LXF A677 tumour growth model parameters.
#


import numpy as np
import pandas as pd
import pints


# Create inverse problem for each mouse ID
mouse_ids = lxf_data['#ID'].unique()
log_likelihoods = []
for ids in mouse_ids:
    # Create mask for mouse with specfied ID
    mouse_mask = lxf_data['#ID'] == ids

    # Get relevant time points
    times = lxf_data[mouse_mask]['TIME in day'].to_numpy()

    # Get measured tumour volumes
    observed_volumes = lxf_data[mouse_mask]['TUMOUR VOLUME in cm^3'].to_numpy()

    # Create inverse problem
    problem = pints.SingleOutputProblem(PintsModel(), times, observed_volumes)

    # Create Gaussian log-likelihood TODO: Change to combined error model
    log_likelihoods.append(pints.GaussianLogLikelihood(problem))

# Create one log_likelihood for the inference from the individual problems
lxf_log_likelihood = pints.SumOfIndependentLogPDFs(log_likelihoods)
#
# Construct likelihood for VXF A341 tumour growth model parameters.
#


import numpy as np
import pandas as pd
import pints


# Create inverse problem for each mouse ID
mouse_ids = vxf_data['#ID'].unique()
log_likelihoods = []
for ids in mouse_ids:
    # Create mask for mouse with specfied ID
    mouse_mask = vxf_data['#ID'] == ids

    # Get relevant time points
    times = vxf_data[mouse_mask]['TIME in day'].to_numpy()

    # Get measured tumour volumes
    observed_volumes = vxf_data[mouse_mask]['TUMOUR VOLUME in cm^3'].to_numpy()

    # Create inverse problem
    problem = pints.SingleOutputProblem(PintsModel(), times, observed_volumes)

    # Create Gaussian log-likelihood TODO: Change to combined error model
    log_likelihoods.append(pints.GaussianLogLikelihood(problem))

# Create one log_likelihood for the inference from the individual problems
vxf_log_likelihood = pints.SumOfIndependentLogPDFs(log_likelihoods)
#
# Check identifiability of LXF A677 problem.
#

import numpy as np
import pints


# Get number of optimisation runs and number of parameters
n = 5
n_parameters = lxf_log_likelihood.n_parameters()

# Initial guess of parameters [initial volume, lambda_0, lambda_1, sigma]
log_parameters_and_sigma = np.array([0.0, 0.0, 0.0, 1.0])

# Standard deviatoin around initial guesses
sigma0 = np.array([0.5, 0.5, 0.5, 0.5])

# Create estimate container
estimates = np.empty(shape=(n, n_parameters))

print('Running...')
# estimate parameters
for i in range(n):
    print('Round %d' % (i+1))

    # Choose optimisation method
    optimiser = pints.CMAES

    # Create optimisation object
    opt = pints.OptimisationController(
        function=lxf_log_likelihood,
        x0=log_parameters_and_sigma,
        sigma0=sigma0,
        method=optimiser)

    # Disable logging mode
    opt.set_log_to_screen(False)

    # Parallelise optimisation
    opt.set_parallel(True)

    # Optimise likelihood
    log_est, _ = opt.run()

    # Save estimates (transform log-parameters, but keep sigma)
    estimates[i, :-1] = np.exp(log_est[:-1])
    estimates[i, -1] =log_est[-1]

print('Done!')

print(' ')
print('Estimates: ')
print('Initial tumour volume [cm^3]: ', estimates[:, 0])
print('Exponential growth rate \lambda _0 [1/day]: ', estimates[:, 1])
print('Linear growth rate \lambda _1 [cm^3/day]: ', estimates[:, 2])
print('Standard deviation of base-level noise [cm^3]: ', estimates[:, 3])<jupyter_output>Running...
Round 1
Round 2
Round 3
Round 4
Round 5
Done!
 
Estimates: 
Initial tumour volume [cm^3]:  [1.34771772 1.34795312 1.34770143 1.34740049 1.34745836]
Exponential growth rate \lambda _0 [1/day]:  [2.03554194e+03 3.66076628e+02 4.82434343e+13 8.42074049e+02
 8.29368456e+02]
Linear growth rate \lambda _1 [cm^3/day]:  [0.35874319 0.35878604 0.358735   0.35877458 0.35877169]
Standard-deviation of base-level noise [cm^3]:  [4.51404352 4.51403982 4.51404453 4.51404477 4.51404398]<jupyter_text>Non-dimensionalise parameters<jupyter_code>#
# Define pints model wrapper such that myokit model can be used for inference.
#
# Input parameters are scaled to an expected scale and log-trasnformed to stabilise
# inference.
#

import myokit
import numpy as np
import pints

from pkpd import model as model


# Wrap myokit model, so it can be used with pints
class PintsModel(pints.ForwardModel):
    def __init__(self):
        # Create myokit model
        model = m.create_tumour_growth_model()

        # Create simulator
        self.sim = myokit.Simulation(model)

        # Characteristic scale [intial volume, lambda_0, lambda_1]
        self._char_scale = np.array([1.3, 3E02, 0.35])

    def n_parameters(self):
        """
        Number of parameters to fit. Here initial V^s_T, lambda_0, lambda_1
        """
        return 3

    def n_outputs(self):
        return 1

    def simulate(self, log_parameters, times):
        # Reset simulation
        self.sim.reset()

        # Sort input parameters
        initial_volume, lambda_0, lambda_1 = np.exp(log_parameters) * self._char_scale

        # Set initial condition
        self.sim.set_state([initial_volume])

        # Set growth constants
        self.sim.set_constant('central.lambda_0', lambda_0)
        self.sim.set_constant('central.lambda_1', lambda_1)

        # Define logged variable
        loggedVariable = 'central.volume_t'

        # Simulate
        output = self.sim.run(times[-1] + 1, log=[loggedVariable], log_times=times)
        result = output[loggedVariable]

        return np.array(result)
#
# Construct likelihood for LXF A677 tumour growth model parameters.
#


import numpy as np
import pandas as pd
import pints


# Create inverse problem for each mouse ID
mouse_ids = lxf_data['#ID'].unique()
log_likelihoods = []
for ids in mouse_ids:
    # Create mask for mouse with specfied ID
    mouse_mask = lxf_data['#ID'] == ids

    # Get relevant time points
    times = lxf_data[mouse_mask]['TIME in day'].to_numpy()

    # Get measured tumour volumes
    observed_volumes = lxf_data[mouse_mask]['TUMOUR VOLUME in cm^3'].to_numpy()

    # Create inverse problem
    problem = pints.SingleOutputProblem(PintsModel(), times, observed_volumes)

    # Create Gaussian log-likelihood TODO: Change to combined error model
    log_likelihoods.append(pints.GaussianLogLikelihood(problem))

# Create one log_likelihood for the inference from the individual problems
lxf_log_likelihood = pints.SumOfIndependentLogPDFs(log_likelihoods)#
# Construct likelihood for LXF A677 tumour growth model parameters.
#


import numpy as np
import pandas as pd
import pints


# Create inverse problem for each mouse ID
mouse_ids = lxf_data['#ID'].unique()
log_likelihoods = []
for ids in mouse_ids:
    # Create mask for mouse with specfied ID
    mouse_mask = lxf_data['#ID'] == ids

    # Get relevant time points
    times = lxf_data[mouse_mask]['TIME in day'].to_numpy()

    # Get measured tumour volumes
    observed_volumes = lxf_data[mouse_mask]['TUMOUR VOLUME in cm^3'].to_numpy()

    # Create inverse problem
    problem = pints.SingleOutputProblem(PintsModel(), times, observed_volumes)

    # Create Gaussian log-likelihood TODO: Change to combined error model
    log_likelihoods.append(pints.GaussianLogLikelihood(problem))

# Create one log_likelihood for the inference from the individual problems
lxf_log_likelihood = pints.SumOfIndependentLogPDFs(log_likelihoods)
#
# Construct likelihood for VXF A341 tumour growth model parameters.
#


import numpy as np
import pandas as pd
import pints


# Create inverse problem for each mouse ID
mouse_ids = vxf_data['#ID'].unique()
log_likelihoods = []
for ids in mouse_ids:
    # Create mask for mouse with specfied ID
    mouse_mask = vxf_data['#ID'] == ids

    # Get relevant time points
    times = vxf_data[mouse_mask]['TIME in day'].to_numpy()

    # Get measured tumour volumes
    observed_volumes = vxf_data[mouse_mask]['TUMOUR VOLUME in cm^3'].to_numpy()

    # Create inverse problem
    problem = pints.SingleOutputProblem(PintsModel(), times, observed_volumes)

    # Create Gaussian log-likelihood TODO: Change to combined error model
    log_likelihoods.append(pints.GaussianLogLikelihood(problem))

# Create one log_likelihood for the inference from the individual problems
vxf_log_likelihood = pints.SumOfIndependentLogPDFs(log_likelihoods)
#
# Check identifiability of LXF A677 problem.
#

import numpy as np
import pints


# Get number of optimisation runs and number of parameters
n = 5
n_parameters = lxf_log_likelihood.n_parameters()

# Initial guess of parameters [initial volume, lambda_0, lambda_1, sigma]
log_parameters_and_sigma = np.array([0.0, 0.0, 0.0, 1.0])

# Standard deviatoin around initial guesses
sigma0 = np.array([0.5, 0.5, 0.5, 0.5])

# Create estimate container
estimates = np.empty(shape=(n, n_parameters))

print('Running...')
# estimate parameters
for i in range(n):
    print('Round %d' % (i+1))

    # Choose optimisation method
    optimiser = pints.CMAES

    # Create optimisation object
    opt = pints.OptimisationController(
        function=lxf_log_likelihood,
        x0=log_parameters_and_sigma,
        sigma0=sigma0,
        method=optimiser)

    # Disable logging mode
    opt.set_log_to_screen(False)

    # Parallelise optimisation
    opt.set_parallel(True)

    # Optimise likelihood
    est, _ = opt.run()

    # Transform parameters back and save estimates
    char_scale = np.array([1.3, 3E02, 0.35])  # Defined in PintsModel
    estimates[i, :-1] = np.exp(est[:-1]) * char_scale
    estimates[i, -1] = est[-1]

print('Done!')

print(' ')
print('Estimates: ')
print('Initial tumour volume [cm^3]: ', estimates[:, 0])
print('Exponential growth rate \lambda _0 [1/day]: ', estimates[:, 1])
print('Linear growth rate \lambda _1 [cm^3/day]: ', estimates[:, 2])
print('Standard deviation of base-level noise [cm^3]: ', estimates[:, 3])<jupyter_output>Running...
Round 1
Round 2
Round 3
Round 4
Round 5
Done!
 
Estimates: 
Initial tumour volume [cm^3]:  [1.34797863 1.34796979 1.34746332 1.3475335  1.34739554]
Exponential growth rate \lambda _0 [1/day]:  [366.0679323  366.07135534 829.36515154 182.86743461 842.07710061]
Linear growth rate \lambda _1 [cm^3/day]:  [0.35878428 0.35878532 0.35877156 0.35883931 0.35877458]
Standard deviation of base-level noise [cm^3]:  [4.5140468  4.51403935 4.51404432 4.51405563 4.51404329]
{'fraction_non_alphanumeric': 0.0969461559067774, 'fraction_numerical': 0.028422180551834986, 'mean_word_length': 3.4584975516541263, 'pattern_counts': {'":': 8, '<': 74, '<?xml version=': 0, '>': 74, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '578064', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 12465, 'n_tokens_neox': 11424, 'n_words': 3685, 'path': 'notebooks/control_growth/identifiability.ipynb', 'repo_name': 'DavAug/ErlotinibGefinitib', 'size': '37330'}