Skip to contents

Systems Biology Models are usually described with reactions, kinetic laws, parameters, and possibly inputs and outputs.

In this R package, we use plain TSV files to describe these quantities. These TSV files have very few rules. One universal rule is that the first column is always reserved for row-names. Row-names must work as variable names in R and C. Certan names would conflict with builtin commands (if, while, do, continue, break). Obviously such names cannot be used for reacting compounds or parameters. In our examples, we label the first column with id (but it can have any name).

They can be imported into R using the command uqsa::model_from_tsv():

library(uqsa)
library(errors)

## Read the TSV files, create an ODE model
m <- model_from_tsv()    # defaults to files from $(pwd)
o <- as_ode(m)           # the ODE

C <- generateCode(o)     # list of strings
odeModel <- write_and_compile(C)
ex <- experiments(m,o)

The model tables follow these guidelines:

  • A list of quantities, like parameters have these columns:
    • id (first column, any heading)
    • value
    • unit
    • scale (optional, assumed to be linear)
  • A list of mathematical expressions can use the term formula instead of value:
    • id
    • formula
    • unit
  • Any quantity that is defined using a row in some table can get updated values in another table, or values to comapre against
  • Data tables need to specifiy values using a notation that keeps values and their uncertainties together, e.g.:
    • 1.234(21)
    • 1.5(1)e+4
  • The tables are found by name, the names our import scripts will be looking for are:
    • Constant.tsv
    • Parameter.tsv
    • Input.tsv
    • Expression.tsv
    • Reaction.tsv
    • Output.tsv
    • Experiment.tsv
  • The Experiment table describes the experimental conditions the model needs to be subjected to in order to replicate measured data:
    • each experimental setting gets a row
    • each row id is used to find the data
    • the data is found in the TSV file named like the row-id
    • Data for experiment EX1 is expected to be in the file EX1.tsv

This package will assume that the directory/folder where the TSV files reside in is the name of the model. This is an example folder structure:

~/Documents/AKAP79$ ls -1 *.tsv *.md
cAMP_increase.tsv
Compound.tsv
EX00___cAMP___CaN_AKAP79.md
EX11____0nM__TRUE___TRUE.tsv
EX12____0nM__TRUE__FALSE.tsv
EX13____0nM_FALSE__FALSE.tsv
EX21__100nM__TRUE___TRUE.tsv
EX22__100nM__TRUE__FALSE.tsv
EX23__100nM_FALSE__FALSE.tsv
EX31__200nM__TRUE___TRUE.tsv
EX32__200nM__TRUE__FALSE.tsv
EX33__200nM_FALSE__FALSE.tsv
EX41__500nM__TRUE___TRUE.tsv
EX42__500nM__TRUE__FALSE.tsv
EX43__500nM_FALSE__FALSE.tsv
EX51_1000nM__TRUE___TRUE.tsv
EX52_1000nM__TRUE__FALSE.tsv
EX53_1000nM_FALSE__FALSE.tsv
EX61_2000nM__TRUE___TRUE.tsv
EX62_2000nM__TRUE__FALSE.tsv
EX63_2000nM_FALSE__FALSE.tsv
EX71___dose__TRUE___TRUE.tsv
EX72__event__TRUE___TRUE.tsv
Experiments.tsv
Expression.tsv
Input.tsv
Output.tsv
Parameter.tsv
Reaction.tsv
README.md
Transformation.tsv

Here is an example command to view a TSV file nicely formatted:

~/Documents/AKAP79$ column -t -s $'\t' Reaction.tsv
id   kinetic law                                reactants        arw  products        is reversible
r51  k5_1*Rii_C                                 Rii_C            <=>  RiiP_C          0
r14  k4_1*RiiP*C - k1_4*RiiP_C                  RiiP + C         <=>  RiiP_C          1
r12  k1_2*RiiP_C*cAMP - k2_1*RiiP_C_cAMP        RiiP_C + cAMP    <=>  RiiP_C_cAMP     1
r43  k4_3*cAMP*RiiP - k3_4*RiiP_cAMP            cAMP + RiiP      <=>  RiiP_cAMP       1
r23  k3_2*RiiP_cAMP*C - k2_3*RiiP_C_cAMP        RiiP_cAMP + C    <=>  RiiP_C_cAMP     1
r78  k8_7*cAMP*Rii - k7_8*Rii_cAMP              cAMP + Rii       <=>  Rii_cAMP        1
r56  k5_6*Rii_C*cAMP - k6_5*Rii_C_cAMP          Rii_C + cAMP     <=>  Rii_C_cAMP      1
r76  k7_6*Rii_cAMP*C - k6_7*Rii_C_cAMP          Rii_cAMP + C     <=>  Rii_C_cAMP      1
r62  k6_2*Rii_C_cAMP                            Rii_C_cAMP       <=>  RiiP_C_cAMP     0
r58  k8_5*Rii*C - k5_8*Rii_C                    Rii + C          <=>  Rii_C           1
r44  k4_4p*RiiP*CaN - k4p_4*RiiP_CaN            RiiP + CaN       <=>  RiiP_CaN        1
r33  k3_3p*CaN*RiiP_cAMP - k3p_3*RiiP_cAMP_CaN  CaN + RiiP_cAMP  <=>  RiiP_cAMP_CaN   1
r48  k4p_8*RiiP_CaN                             RiiP_CaN         <=>  Rii + CaN       0
r37  k3p_7*RiiP_cAMP_CaN                        RiiP_cAMP_CaN    <=>  CaN + Rii_cAMP  0
r1   kf_C_AKAR4*C*AKAR4 - kb_C_AKAR4*AKAR4_C    C + AKAR4        <=>  AKAR4_C         1
r2   kcat_AKARp*AKAR4_C                         AKAR4_C          <=>  AKAR4p + C      0

In the following sections, we give examples for each of the tables, as markdown tables.

Reaction

id kinetic law reactants arw products is reversible
r51 k5_1*Rii_C Rii_C <=> RiiP_C 0
r14 k4_1*RiiP*C - k1_4*RiiP_C RiiP + C <=> RiiP_C 1
r12 k1_2*RiiP_C*cAMP - k2_1*RiiP_C_cAMP RiiP_C + cAMP <=> RiiP_C_cAMP 1
r43 k4_3*cAMP*RiiP - k3_4*RiiP_cAMP cAMP + RiiP <=> RiiP_cAMP 1
r23 k3_2*RiiP_cAMP*C - k2_3*RiiP_C_cAMP RiiP_cAMP + C <=> RiiP_C_cAMP 1
r78 k8_7*cAMP*Rii - k7_8*Rii_cAMP cAMP + Rii <=> Rii_cAMP 1
r56 k5_6*Rii_C*cAMP - k6_5*Rii_C_cAMP Rii_C + cAMP <=> Rii_C_cAMP 1
r76 k7_6*Rii_cAMP*C - k6_7*Rii_C_cAMP Rii_cAMP + C <=> Rii_C_cAMP 1
r62 k6_2*Rii_C_cAMP Rii_C_cAMP <=> RiiP_C_cAMP 0
r58 k8_5*Rii*C - k5_8*Rii_C Rii + C <=> Rii_C 1
r44 k4_4p*RiiP*CaN - k4p_4*RiiP_CaN RiiP + CaN <=> RiiP_CaN 1
r33 k3_3p*CaN*RiiP_cAMP-k3p_3*RiiP_cAMP_CaN CaN + RiiP_cAMP <=> RiiP_cAMP_CaN 1
r48 k4p_8*RiiP_CaN RiiP_CaN <=> Rii + CaN 0
r37 k3p_7*RiiP_cAMP_CaN RiiP_cAMP_CaN <=> CaN + Rii_cAMP 0
r1 kf_C_AKAR4*C*AKAR4-kb_C_AKAR4*AKAR4_C C + AKAR4 <=> AKAR4_C 1
r2 kcat_AKARp*AKAR4_C AKAR4_C <=> AKAR4p + C 0

Compound

Reacting compounds, or molecular species. These will be turned into state variables.

id value Unit uniprot Comment
Rii 6.3 µM KAP2_HUMAN also KAP3_HUMAN because Rii = RII_alpha+RII_beta
cAMP 0 µM NONE NONE
RiiP 0 µM NONE phosphorylated
Rii_C 0.63 µM NONE NONE
RiiP_cAMP 0 µM NONE NONE
RiiP_C 0 µM NONE NONE
RiiP_C_cAMP 0 µM NONE NONE
C 0 µM KAPCB_HUMAN NONE
Rii_cAMP 0 µM NONE NONE
Rii_C_cAMP 0 µM NONE NONE
CaN 1.5 µM PP2BA_HUMAN PP3R1_HUMAN
RiiP_CaN 0 µM NONE NONE
RiiP_cAMP_CaN 0 µM NONE NONE
AKAR4 0.2 µM NONE RRID:addgene_61619
AKAR4_C 0 µM NONE NONE
AKAR4p 0 µM NONE NONE

Parameter

These are parameters, that will be subject to optimization or sampling. The values are assumed to be crude estimates.

Known parameters should go into a different table.

id value median stdv prior distribution scale unit
k5_1 1.667837 1.519 3 Normal distribution log10 s^-1
k1_2 0.0946379 −0.306 3 Normal distribution log10 1/µM*s
k3_2 −2.483545 −2.264 3 Normal distribution log10 1/µM*s
k2_3 −0.230731 −1.807 3 Normal distribution log10 s^-1
k3_4 −1.136198 −2.796 3 Normal distribution log10 s^-1
k4_3 −2.147382 −1.824 3 Normal distribution log10 1/µM*s
k4_1 −2.284185 −1.42 3 Normal distribution log10 1/µM*s
k1_4 −2.928144 −2.585 3 Normal distribution log10 s^-1
k8_7 −2.149449 −1.824 3 Normal distribution log10 1/µM*s
k7_8 −1.07931 −2.796 3 Normal distribution log10 s^-1
k5_6 −0.177019 −0.305 3 Normal distribution log10 1/µM*s
k6_5 −0.225915 0.15 3 Normal distribution log10 s^-1
k8_5 −0.998736 0.322 3 Normal distribution log10 1/µM*s
k7_6 −1.00916 −0.525 3 Normal distribution log10 1/µM*s
k6_7 −1.653377 −1.745 3 Normal distribution log10 s^-1
k6_2 1.5293947 1.519 3 Normal distribution log10 s^-1
k5_8 −3.696929 −3.523 3 Normal distribution log10 s^-1
AKAPoff_1 −0.531584 0.415 3 Normal distribution log10 s^-1
AKAPoff_3 1.1736716 1.301 3 Normal distribution log10 s^-1
AKAPon_1 −0.365851 −0.347 0.279 Normal distribution log10 s^-1
AKAPon_3 0.2361535 0.301 3 Normal distribution log10 s^-1
kf_C_AKAR4 −1.74248 −1.745 0.0969 Normal distribution log10 1/µM*s
kb_C_AKAR4 −0.98046 −0.975 0.0969 Normal distribution log10 s^-1
kcat_AKARp 1.0077982 1.009 0.0969 Normal distribution log10 s^-1
kmOFF 2.0096025 2 0.176 Normal distribution log10 µM
kmON −0.005704 0 0.176 Normal distribution log10 µM
KD_T −0.175253 −0.155 0.301 Normal distribution log10 µM

Input

These are known parameter values. Parameters that are set by us in exactly the same way the experimentalost controls the conditions of the experiment.

id value unit
b_AKAP 0 1

Expression

These are mathematical expression that can be used as abbreviations, re-used quantities that perhaps appear in several kinetic laws.

id formula unit
k3p_7 b_AKAP*AKAPon_1 + (1 - b_AKAP)*AKAPoff_1 1/s
k4p_4 b_AKAP*AKAPon_3 + (1 - b_AKAP)*AKAPoff_3 1/s
k4p_8 b_AKAP*AKAPon_1 + (1 - b_AKAP)*AKAPoff_1 1/s
k3p_3 b_AKAP*AKAPon_3 + (1 - b_AKAP)*AKAPoff_3 1/s
k4_4p b_AKAP*(k4p_4 + k3p_7)/kmON + (1 - b_AKAP)*(k4p_4 + k3p_7)/kmOFF 1/(micromolarity*s)
k3_3p b_AKAP*(k3p_3 + k4p_8)/kmON + (1 - b_AKAP)*(k3p_3 + k4p_8)/kmOFF 1/(micromolarity*s)
k2_1 k1_2 * KD_T 1/s

Experiment

This is the table of experiments. Each column can name an input or compound, these value smodify the initial conditions (and inputs) for that experiment.

id cAMP CaN b_AKAP t0 type event comment
EX11____0nM__TRUE___TRUE 0 1.5 1 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX12____0nM__TRUE__FALSE 0 1.5 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX13____0nM_FALSE__FALSE 0 0 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX21__100nM__TRUE___TRUE 0.1 1.5 1 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX22__100nM__TRUE__FALSE 0.1 1.5 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX23__100nM_FALSE__FALSE 0.1 0 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX31__200nM__TRUE___TRUE 0.2 1.5 1 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX32__200nM__TRUE__FALSE 0.2 1.5 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX33__200nM_FALSE__FALSE 0.2 0 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX41__500nM__TRUE___TRUE 0.5 1.5 1 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX42__500nM__TRUE__FALSE 0.5 1.5 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX43__500nM_FALSE__FALSE 0.5 0 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX51_1000nM__TRUE___TRUE 1 1.5 1 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX52_1000nM__TRUE__FALSE 1 1.5 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX53_1000nM_FALSE__FALSE 1 0 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX61_2000nM__TRUE___TRUE 2 1.5 1 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX62_2000nM__TRUE__FALSE 2 1.5 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX63_2000nM_FALSE__FALSE 2 0 0 −30 time series Figure 3 in https://doi.org/10.7554/eLife.68164
EX71___dose__TRUE___TRUE 0 1.5 1 −30 dose response none
EX72__event__TRUE___TRUE 0 1.5 1 −30 time series cAMP_increase none

Output

Each row describes an observable, and how to measure it with this model (as a function of state variables).

id formula unit
AKAR4pOUT (AKAR4p*5)*71.67+100 µM

Transformation

Each row describes a transformation, each column names the quantity to be transformed, any input and compound can be named here. Each cell contains a formula that calculates the new value. The id can be used in event schedules to reference this transformation.

id cAMP
SET_cAMP dose

Event

An experiment can have an event schedule: a timed list of sudden transformation, each event occurs at a given time, applies a transformation (by name) and assigns a _dose_, the intensity of the transformation.

id time transformation dose
tr1 610 SET_cAMP 0.1
tr2 810 SET_cAMP 0.2
tr3 1010 SET_cAMP 0.5
tr4 1210 SET_cAMP 1.0
tr5 1410 SET_cAMP 2.0

Data

This is an example data-set. Each measured value is given with the uncertainty in the same cell. The value 9.994(45)E1 means 99.94 with standard error of 0.45. Sometimes this is written as 99.94±0.45 (whichthe package also accepts). But it’s the inferior syntax for this as it collides with specifying two solutions of some sort, e.g.: 1 ± sqrt(2).

A Time Series Experiment

id time AKAR4pOUT
E0103T001 −15 9.994(45)E1
E0103T002 −10 9.997(88)E1
E0103T003 −5 9.960(31)E1
E0103T004 0 NA
E0103T005 5 9.80(14)E1
E0103T006 10 9.81(19)E1
E0103T007 15 9.701(98)E1
E0103T008 20 9.79(15)E1
E0103T009 25 9.86(11)E1
E0103T010 30 9.95(14)E1
E0103T011 35 9.87(15)E1
E0103T012 40 9.91(14)E1
E0103T013 45 9.93(13)E1
E0103T014 50 9.87(13)E1
E0103T015 55 9.815(68)E1
E0103T016 60 9.8(1)E1
E0103T017 65 9.94(11)E1
E0103T018 70 9.73(11)E1
E0103T019 75 9.769(43)E1
E0103T020 80 9.95(9)E1
E0103T021 85 9.90(14)E1
E0103T022 90 9.7(1)E1
E0103T023 95 9.833(88)E1
E0103T024 100 9.812(29)E1
E0103T025 105 9.883(77)E1
E0103T026 110 1.006(11)E2
E0103T027 115 9.960(78)E1
E0103T028 120 9.91(12)E1
E0103T029 125 1.006(14)E2
E0103T030 130 9.81(15)E1
E0103T031 135 9.99(14)E1
E0103T032 140 9.9(1)E1
E0103T033 145 9.87(15)E1
E0103T034 150 9.80(16)E1
E0103T035 155 9.812(85)E1
E0103T036 160 9.9(1)E1
E0103T037 165 9.948(88)E1
E0103T038 170 9.85(12)E1
E0103T039 175 9.88(14)E1
E0103T040 180 9.90(13)E1
E0103T041 185 9.85(22)E1
E0103T042 190 9.963(71)E1
E0103T043 195 9.948(94)E1
E0103T044 200 9.858(65)E1
E0103T045 205 9.966(69)E1
E0103T046 210 9.88(8)E1
E0103T047 215 9.969(59)E1
E0103T048 220 9.972(14)E1
E0103T049 225 9.809(87)E1
E0103T050 230 9.90(13)E1
E0103T051 235 9.89(14)E1
E0103T052 240 9.96(12)E1
E0103T053 245 9.79(15)E1
E0103T054 250 9.7(2)E1
E0103T055 255 9.849(48)E1
E0103T056 260 9.95(24)E1
E0103T057 265 9.935(86)E1
E0103T058 270 9.82(16)E1
E0103T059 275 9.94(4)E1
E0103T060 280 9.84(13)E1
E0103T061 285 9.97(14)E1
E0103T062 290 9.966(45)E1
E0103T063 295 9.79(14)E1
E0103T064 300 9.93(12)E1
E0103T065 305 9.843(14)E1
E0103T066 310 9.97(17)E1
E0103T067 315 9.84(13)E1
E0103T068 320 9.8(1)E1
E0103T069 325 9.88(12)E1
E0103T070 330 9.827(47)E1
E0103T071 335 9.83(12)E1
E0103T072 340 9.78(16)E1
E0103T073 345 9.988(86)E1
E0103T074 350 10.0(1)E1
E0103T075 355 9.972(51)E1
E0103T076 360 9.923(57)E1
E0103T077 365 9.96(11)E1
E0103T078 370 9.9(1)E1
E0103T079 375 9.88(15)E1
E0103T080 380 9.9(1)E1
E0103T081 385 9.806(28)E1
E0103T082 390 1.0000(61)E2
E0103T083 395 9.92(15)E1
E0103T084 400 9.91(14)E1
E0103T085 405 1.009(14)E2
E0103T086 410 1.0019(77)E2
E0103T087 415 9.985(99)E1
E0103T088 420 1.0000(92)E2
E0103T089 425 9.86(15)E1
E0103T090 430 9.877(48)E1
E0103T091 435 9.99(14)E1
E0103T092 440 9.892(44)E1
E0103T093 445 9.92(13)E1
E0103T094 450 1.006(14)E2
E0103T095 455 1.0086(98)E2
E0103T096 460 9.96(12)E1
E0103T097 465 9.95(11)E1
E0103T098 470 9.932(99)E1
E0103T099 475 9.99(12)E1
E0103T100 480 1.0049(22)E2
E0103T101 485 9.9(1)E1
E0103T102 490 10.0(1)E1
E0103T103 495 1.00(1)E2
E0103T104 500 9.98(12)E1
E0103T105 505 1.0114(32)E2
E0103T106 510 1.007(12)E2
E0103T107 515 9.972(45)E1
E0103T108 520 1.0025(92)E2
E0103T109 525 1.010(15)E2
E0103T110 530 9.98(2)E1
E0103T111 535 1.0114(79)E2
E0103T112 540 1.004(16)E2
E0103T113 545 9.96(4)E1
E0103T114 550 1.0034(82)E2
E0103T115 555 9.95(18)E1
E0103T116 560 1.0022(68)E2
E0103T117 565 9.89(24)E1
E0103T118 570 9.963(69)E1
E0103T119 575 1.010(11)E2
E0103T120 580 9.94(12)E1
E0103T121 585 9.938(92)E1
E0103T122 590 1.0090(94)E2
E0103T123 595 9.94(14)E1
E0103T124 600 1.008(14)E2
E0103T125 605 1.0028(73)E2

There are many tables like this, one per experiment.

A Dose Response Experiment

A Dose Response table is a more convenient way of describing a series of experiments, each with its own, very _similar_ conditions, and exactly one measurement per time-course (a series of micro experiments). Such a table will be split up by our import script, and one additional simulation experiment per line will be created. The newly created simulation-experiments will have very similar names and vcan be found this way.

id cAMP time AKAR4pOUT
E0103T125 0 605 1.0028(73)E2
E0203T125 0.1 605 1.0003(91)E2
E0403T125 0.2 605 1.10(6)E2
E0503T125 0.5 605 1.33(12)E2
E0603T125 1 605 1.5336(97)E2
E0303T125 2 605 1.675(37)E2

If you don’t want these values to be split up into several micro-experiments, then a very similar effect can be achieved using events and transformations. See below.

An Event Schedule Experiment

This is one time series experiment that replicates the same data-set as the dose reponse curve in the previous section. But, in this sequence, we also have sudden events that change the cAMP value in the system, and then we allow the model to find a new equilibrium (in 200 seconds).

id time AKAR4pOUT
EVENTEXT01 605 1.0028(73)E2
EVENTEXT02 805 1.0003(91)E2
EVENTEXT03 1005 1.10(6)E2
EVENTEXT04 1205 1.33(12)E2
EVENTEXT05 1405 1.5336(97)E2
EVENTEXT06 1605 1.675(37)E2

The event schedule we used here is the same as in the section about events:

id time transformation dose
tr1 610 SET_cAMP 0.1
tr2 810 SET_cAMP 0.2
tr3 1010 SET_cAMP 0.5
tr4 1210 SET_cAMP 1.0
tr5 1410 SET_cAMP 2.0

These two work together to create a sequence of letting the moment converge to an equilibrium, recording the value of AKAR4pOUT, then setting a new cAMP level and repeat.