model_building.Rmd

---
title: "Model Building"
author: "Jonas Förster & Jürgen Pahle"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Model Building}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, echo=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.asp = 10/16,
  out.width = "100%",
  error = TRUE,
  warning = TRUE,
  message = FALSE
)
set.seed(1)
```

# Initial Setup

```{r setup}
library(ggplot2)
library(CoRC)
```

# Quick Model Building

Building a quick Michaelis-Menten reaction using COPASI's automatisms.
Working with reaction equations automatically creates required entities the same way the COPASI GUI does.

```{r}
newModel()

newReaction("a = b", fun = "Reversible Michaelis-Menten")

getSpecies()
getReactions()
getParameters()

unloadModel()
```

# More Detail

Getting more detail in can be helpful.
Here we create a global quantity "substrate concentration" which applies a fixed value to species "a".
This is achieved by assigning the `Value` reference of the global quantity to species "a".

```{r}
newModel()

subs_quant <- newGlobalQuantity("substrate concentration", initial_value = 10)
subs_value <- quantity(subs_quant, reference = "Value")
subs_value

newSpecies("a", type = "assignment", expression = subs_value)
newReaction("a = b", fun = "Reversible Michaelis-Menten")

getGlobalQuantities()
getSpecies()
getReactions()
getParameters()

timecourse <- runTimeCourse(1000)

autoplot.copasi_ts(timecourse)

unloadModel()
```

# Defining Reaction Kinetics Manually

In previous examples, we used the predefined "Reversible Michaelis-Menten" kinetic functions.
In this example we create an equivalent function manually to achieve the same result.

```{r}
newModel()

newCompartment("compartment")
newSpecies("a", type = "fixed", initial_concentration = 10)
newSpecies("b")
newKineticFunction(
  "manual MM",
  "((Vf * substrate) / Kms - Vr * product / Kmp) / (1 + substrate / Kms + product / Kms)",
  parameters = list(substrate = "substrate", product = "product")
)
newReaction(
  "a = b",
  fun = "manual MM",
  mapping = list(
    substrate = "a",
    product = "b",
    Vf = 0.1,
    Kms = 0.1,
    Vr = 0.1,
    Kmp = 0.1
  )
)

timecourse <- runTimeCourse(1000)

autoplot.copasi_ts(timecourse)

unloadModel()
```

# Quantities as Proxy for Time Course Output

The `runTimeCourse` function only returns data for actual defined entities of the model.
Auxiliarry information like reaction rates can be missing and needs to be added via global quantity assignments.
This example creates a global quantity "reaction rate" to generate data for the main reaction's flux.

```{r}
newModel()

subs_quant <- newGlobalQuantity("substrate concentration", initial_value = 10)
subs_value <- quantity(subs_quant, reference = "Value")

newSpecies("a", type = "assignment", expression = subs_value)
newReaction("a = b", fun = "Reversible Michaelis-Menten")

runTC(1000)$result

newGlobalQuantity(
  "reaction rate",
  type = "assignment",
  expression = reaction("a = b", reference = "Flux")
)

runTimeCourse(1000)$result

unloadModel()
```