# Simulates Network Dynamics under Iterative Perturbations 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. ```r perturbation.thread( system, M, x0, perturbation_function, initial_parameter_function, times = seq(0, 100, 1), reduction = identity, to.numeric = FALSE, only.final.state = TRUE ) ``` ## 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. ## Returns 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. ## Examples ```r 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) ```

perturbation.thread function