Rato (version 0.1.0)

perturbation.thread: Simulates Network Dynamics under Iterative Perturbations

Description

This function simulates the behavior of a network undergoing multiple iterations of custom perturbations. Starting with an initial healthy network, it calculates the network's trajectory by solving the system of ordinary differential equations (ODEs) after each perturbation. Perturbations are applied iteratively, followed by dimension reduction using the specified reduction function.

Usage

perturbation.thread(
  system,
  M,
  x0,
  perturbation_function,
  initial_parameter_function,
  times = seq(0, 100, 1),
  reduction = identity,
  to.numeric = FALSE,
  only.final.state = TRUE
)

Value

Depending on `to.numeric`, returns either a list or a numeric matrix, representing the system's state across multiple perturbations. If `to.numeric` is `FALSE`, the function returns a list `L`, where `L[[i]]` represents the final state of the system after the `i`-th perturbation. If `TRUE`, the list is converted to a numeric matrix before being returned.

Arguments

system

A function defining the system's dynamics: `system(time, x, params)` which returns a list of state deltas `dx`.

M

The initial adjacency matrix of the network.

x0

Initial conditions of the network's nodes (numeric vector).

perturbation_function

A function defining the perturbation applied to the network. Its signature should be `perturbation(params, iter)`, where `params` are the current parameters passed to `system`, and `iter` is the current iteration. Access the current network state via `params$M`. It should return updated parameters or `NULL` if no further perturbations should be applied.

initial_parameter_function

A function of type `f(M) -> list` that takes the adjacency matrix of the network as input and returns the initial parameters of the system.

times

A vector of time points for the ODE solver.

reduction

A reduction function applied to the ODE solution. This can be `identity` (for all node states) or functions like `mean` or `median`. The function signature should be either `f(numeric) -> numeric` or `f(matrix) -> matrix`, depending on whether `only.final.state` is `TRUE` or `FALSE`.

to.numeric

Logical; if `TRUE`, converts the final list to a numeric matrix. If `FALSE`, returns a list `L`, where `L[[i]]` is the network state after `i` perturbations, post-reduction.

only.final.state

Logical; if `TRUE`, applies reduction only to the final state of the system. If `FALSE`, the reduction function is applied to the full ODE trajectory, which can be memory-intensive.

Examples

# \donttest{
   node_file <- system.file("extdata", "IL17.nodes.csv", package = "Rato")
   edge_file <- system.file("extdata", "IL17.edges.csv", package = "Rato")
   g <- Rato::graph.from.csv(node_file, edge_file, sep=",", header=TRUE)
 
   initial_params = list('f' = 1, 'h'=2, 'B'=0.01)
   removal_order = NULL
   update_params = identity
   reduction = identity
   initial_parameter_function <- function(M) {

     if (is.list(initial_params)) {
       params <- initial_params
     }
     else if (is.function(initial_params)) {
       params <- initial_params(M)
     }
     n <- nrow(M) # The number of nodes of the graph

     if(is.null(removal_order)){
         removal_order <- sample(1:n, n) 
     }

     if(is.null(params$M)){
       params$M = M
     } 
     params$removal_order = removal_order
     return(params)
   }
   perturbation <- function(params, iter) {
     removal_order <- params$removal_order
     M <- params$M
 
     if(length(removal_order) > 0 ) {
       index <- removal_order[1] 
       removal_order <- removal_order[-1] 
       M[index, ] <- 0 # Remove entries
       M[, index] <- 0 # Remove entries
 
       params$M = M
       params$removal_order = removal_order
 
       # Add the list of removed indices to the parameter list.
       if(is.null(params$removed_indices)){
         params$removed_indices <- c(index)
       } else {
         removed_indices <- c(params$removed_indices, index)
         params$removed_indices <- removed_indices
       }
 
       output <- update_params(params)
       return(output) # Return updated parameters
     }
 
     # If we removed every single entry, just STOP.
     return(NULL)
   }
 
   Rato::perturbation.thread( Rato::Michaelis.Menten
                             , g$M
                             , g$initial_values
                             , initial_parameter_function = initial_parameter_function
                             , perturbation_function = perturbation)

# }