4. A Priori Identifiability Analysis

In this example, we show you how to a NOT Gate system with BMSS2.

A Priori Identifiability Analysis

We want to know if our candidate models will be fully identifiable based on the states measured. This can achieved using the strike-GOLDD algorithm.

import matplotlib.pyplot as plt
import numpy             as np
import pandas            as pd


import BMSS.models.model_handler    as mh
import BMSS.models.setup_sg         as ssg
import BMSS.models.ia_results       as ir
import BMSS.strike_goldd_simplified as sg

Connected to MBase_models, UBase_models
Connected to MBase_settings, UBase_settings

After importing the required modules, we read the config files to set up the necessarydata structures. Note: Note that in the first example, the settings were put in a separate file(s). This allowed us to switch between different groups of settings. However, since we only need one group of settings in this example, we can combine that group of settings with the core_model file for brevity.

[system_type]
system_type = LogicGate, Not, Single

[states]
states = p1, p2

[parameters]
parameters = kp1, synp1, synp2, rep, degp

[inputs]
inputs = u1

[equations]
equations =
    dp1 = synp1 *u1                    -degp*p1
    dp2 = synp2 *(kp1+rep*p1)/(kp1+p1) -degp*p2



[LogicGate, Not, Single]

init =
    p1 = [0],
    p2 = [0]

guess =
    kp1   = [5e-5],
    synp1 = [1e-5],
    synp2 = [1e-5],
    rep   = [0.8],
    degp  = [0.012],
    u1    = [0]

priors =
    degp = [0.012, 0.04]

input_conditions =
    u1 = 3

fixed_parameters =
    [degp, u1]

measured_states =
    [p1, p2]

decomposition =
    [[p1, p2]]

Decomposition

The symbolic algebraic operations used in STRIKE-GOLDD are complex and large models cannot be analyzed at once. Instead, they need to be decomposed into smaller ones first. The decomposition section contains a list of list of state variables where each inner list is a group of state variables that will be analyzed by the algorithm. Each group must contain at least one measured state variable. The state variables do not need to be split into mutually exclusive groups but can be reused e.g. [[p1], [p1, p2]]. However, to ensure comprehensive coverage, check that each state variable appears at least once in your decomposition.

Input Conditions

The input conditions are the number of different value of inputs minus one. For example, u1 in this model represents inducer concentration. If we characterize the system with 4 inducer concentrations, set u1=3.

#Set up core models and sampler arguments
model_files = ['LogicGate_Not_Single.ini',
               ]

user_core_models = [mh.from_config(filename) for filename in model_files]
user_core_models = {core_model['system_type']: core_model for core_model in user_core_models}

sg_args, config_data, variables = ssg.get_strike_goldd_args(model_files, user_core_models=user_core_models)

sg_results = sg.analyze_sg_args(sg_args)
outfile   = 'sg_results.yaml'
yaml_dict = ir.export_sg_results(sg_results, variables, config_data, user_core_models=user_core_models, filename=outfile)

States :  [p1_1, p2_1]
Outputs:  [p1_1, p2_1]
Params :  [kp1_1, synp1_1, synp2_1, rep_1, degp_1]
Inputs :  {u1_1: 3}
ICS    :  {p1_1: 0.0, p2_1: 0.0}
Equations:
-degp_1*p1_1 + synp1_1*u1_1
-degp_1*p2_1 + synp2_1*(kp1_1 + p1_1*rep_1)/(kp1_1 + p1_1)
Building Oi with 3 derivatives and size (8, 7)
Time for adding one Lie derivative: 0.17090320587158203
Time for adding one Lie derivative: 0.25185346603393555
Time for adding one Lie derivative: 0.32781291007995605
Time for building Oi: 0.8884921073913574
Checking rank for matrix with size (8, 7)
Rank 6
Time for Row Reduction: 10.525646686553955
Elim-recalc for matrix with size (8, 7)
Found [synp1_1, synp2_1, degp_1]
Extending Oi
Time for adding one Lie derivative: 0.5756707191467285
Checking rank for matrix with size (10, 7)
Rank 6
Time for Row Reduction: 1.8969168663024902
Elim-recalc for matrix with size (10, 7)
Found []
14.34646224975586
...................................................
Total Time: 14.348460912704468

The results are returned inside the dictionary dst.

yaml_dict[1]

{'system_type': 'LogicGate, Not, Single',
 'init': {'p1': 0.0, 'p2': 0.0},
 'input': {'p1': 1,
  'p2': 1,
  'kp1': 0,
  'synp1': 0,
  'synp2': 0,
  'rep': 0,
  'degp': 0,
  'u1': 3},
 'output': {'p1': 1,
  'p2': 1,
  'kp1': 0,
  'synp1': 1,
  'synp2': 1,
  'rep': 0,
  'degp': 1}}