Created
July 3, 2024 19:05
-
-
Save lesserfish/c83a30cb1421440aa2889fc1acc4efff to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| library(magrittr) | |
| library(deSolve) | |
| sinput <- function(A, i = NULL) | |
| { | |
| sinput_help <- function(A, i) { | |
| A[i, ] %>% sum() %>% return() | |
| } | |
| if(is.null(i)) | |
| { | |
| n <- nrow(A) | |
| sapply(1:n, function(i) sinput_help(A, i)) %>% return() | |
| } | |
| else | |
| { | |
| sinput_help(A, i) %>% return() | |
| } | |
| } | |
| beff <- function(A) | |
| { | |
| n <- nrow(A) | |
| ones <- rep(1, n) | |
| num <- t(ones) %*% A %*% sinput(A) | |
| den <- t(ones) %*% A %*% ones | |
| if(den == 0){ | |
| den = 1 | |
| } | |
| num / den %>% return() | |
| } | |
| # This is the modified Michaelis-Menten equations according to Barabasi. | |
| dxdt <- function(x, beffA, B, f, h) | |
| { | |
| out <- - B*(x**f) + beffA * ((x**h) / (x**h + 1)) | |
| return(out) | |
| } | |
| demo <- function(t, x, params) | |
| { | |
| beffA <- params$beffA | |
| f <- params$f | |
| h <- params$h | |
| B <- params$B | |
| dxdt(x, beffA, B, f, h) %>% list() %>% return() | |
| } | |
| # Perturbation test | |
| custom_perturbation <- function(A) | |
| { | |
| n <- dim(A)[1] | |
| removal_order <- sample(1:n, replace=FALSE) | |
| perturbation <- function(context) | |
| { | |
| params <- context$params | |
| it <- context$iteration | |
| A <- params$A | |
| outer <- params$outer | |
| # Select which node will be removed | |
| node <- removal_order[it] | |
| # Remove outgoing edges | |
| A[node,] <- 0 | |
| outer[node] <- 1 # If we remove an outer node, it's expression can (and will) drop | |
| # Remove incoming edges | |
| A[, node] <- 0 | |
| params[['A']] <- A | |
| params[['beffA']] <- beff(A) | |
| params[['outer']] <- outer | |
| return(params) | |
| } | |
| return(perturbation) | |
| } | |
| custom_parameters = function(A) | |
| { | |
| parameters <- list('f' = 1, 'h'=2, 'A'=A, 'beffA'=beff(A), 'B'=-1) | |
| return(parameters) | |
| } | |
| load("~/Documents/Work/SysBio2023/KEGG_DATABASE_FULL/ko00360/ko00360.RData") | |
| # Convert it to a matrix | |
| M <- as.matrix(igraph::as_adjacency_matrix(G, names=FALSE)) | |
| data <- perturbation_analysis(M = M, | |
| x0 = 2, | |
| perturbation_function = custom_perturbation, | |
| parameter_function = custom_parameters, | |
| dynamics = demo, | |
| iterations = 100, | |
| log = TRUE) | |
| plot_resilience_drop(data, max_y = 2e18) | |
| # Traditional way | |
| n <- dim(M)[1] | |
| x0 = rep(2, n) | |
| data <- perturbation_analysis(M = M, x0 = x0, iterations = 100, log=TRUE) | |
| plot_resilience_drop(data) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment