Install needed libraries:

library(deSolve)
library(FME)
library(rootSolve)
library(ggplot2)
library(ggplotify)
library(patchwork)
library(modelsummary)
library(ggpubr)
library(tidyr)
library(tidyverse)
library(dplyr)
library("TTR")
library(profileR)
library(gginnards)
library(ggpmisc)
library(directlabels)
library("ggrepel")
library(glue)

Libraries.

For now you need to know that libraries are usefull code, running mathematical and statistical calculations.

There are also libraries for different kind of plotting based on the calculations.

The code

The article: Hoermann, R., Pekker, M. J., Midgley, J. E., Larisch, R., & Dietrich, J. W. (2022). Principles of Endocrine Regulation: Reconciling Tensions Between Robustness in Performance and Adaptation to Change. Frontiers in Endocrinology. https://doi.org/10.3389/fendo.2022.825107 is the base for the free code provided by the authors.

The code for the R-script is here: https://www.frontiersin.org/articles/10.3389/fendo.2022.825107/full#supplementary-material

I have changed a little in way the script is build. Parameters and ranges are placed in the beginning of the script as you need to change these in order to calculate different scenarios.

Parameters

# Make your changes to the system parameters
s1=1  #sys1 s1=1
k134=1  #sys1 k134=1
#THR
k21=1  #sys1 k21=1  <1 decreases TSH
c2=1  #sys1 c2=1
k234=1  #sys1 k234=1
#TSH
k32=1  #sys1 k32=1
#FT4
k43=0.4  #sys1 k43=0, sys10.3 = k43=0.2
k42=0.1  #sys1 k42=0, sys10.3 = k42=0.1
k423=0.7  #sys1 k423=1 sys10.3 k423=0.7
#FT3
drug4=0
drug3=0 
a1=1 
a2=1 
a3=1 
a4=1

As you can see from the comments (#text) different values are noted for different systems tested in the article. In the article there are 12 systems. I have focused on sys1 and sys10.3.

Changes in the parameters are reflected in the calculations and in the plots.

Parameter ranges

# Make your changes to the ranges
#THR
s1min=1 
s1max=1
    k134min=1
    k134max=1
#TSH
k21min=1  
k21max=1
    c2min=1
    c2max=1
k234min=1 
k234max=1
#FT4
k32min=0.3
k32max=3
#FT3
k43min=0.3 #sys10.3 = k43=0.7
k43max=0.4 #sys10.3 = k43=0.8
    k42min=0.3  #sys10.3 = k42=0.1
    k42max=0.3 #sys10.3 = k42=0.2
k423min=0.2 #sys10.3 = k423=0.2 
k423max=0.3 #sys10.3 = k423=0.3
#not necessary as of now
# drug4=0  
# drug3=0 
# a1=1 
# a2=1 
# a3=1 
# a4=1

In the calculations you are able to use ranges for the parameters. This enables you to judge the effect of a broad range of parameters as opposed to evaluating one parameter at a time.

Establishing “k” the combination of parameters used in the calculations.

It is important that you do not change anything in this chunck of code.

# system parameters ----
# DO NOT CHANGE This is only establishing "k"
k=c(s1=s1, #sys1 s1=1
    k134=k134, #sys1 k134=1
    #THR
        k21=k21, #sys1 k21=1, <1 decreases TSH
        c2=c2, #sys1 c2=1
        k234=k234, #sys1 k234=1
        #TSH
            k32=k32, #sys1 k32=1
            #FT4
                k43=k43, #sys1 k43=0
                k42=k42, #sys1 k42=0
                k423=k423, #sys1 k423=1
                #FT3
    drug4=drug4, 
    drug3=drug3,
    a1=a1,
    a2=a2,
    a3=a3,
    a4=a4)

Naming the plots

Give your plots a meaningful title:

#change the plot names
#REMEMBER
p_titles <-c("Sys10.3 - ")

sys10.3 may not be very informative so add your own title.

Show values

Show values different from ONE as a subtitle in the plots.

#Extract those with values different from one
k2 <- k[k != 1]
k2_df = as.data.frame((k2))
#create names and values for k
k2_df$names<-c(names(k2))

Starting the scripts

Scripts to generate the graphical solutions and figures. These lines are from the supplemental material of the article.

fn1=function(t,y,k){ 
    with (as.list(c(y,k)), { 
        dTRH= s1/(k134*FT4*FT3) - a1*TRH 
        dTSH= k21*TRH+c2/(k234*FT4*FT3) - a2*TSH 
        dFT4= k32*TSH + drug4 - a3*FT4 #dFT4= k32*TSH - a3*FT4 ->sys1
        dFT3= k423*TSH*FT4 + k43*FT4 + k42*TSH + drug3 - a4*FT3 #dFT3= k423*TSH*FT4 - a4*FT3 ->sys1
        return(list(c(dTRH,dTSH,dFT4,dFT3))) 
    }) 
} 

# initial values 
y0=c(TRH=1,TSH=1,FT4=1,FT3=1) 
# time-dependent solution 
times=seq(0,20, by=0.1) 
out1=ode(y=y0,func=fn1,times=seq(0,20,by=0.1),parms=k)

Here we create “out1” which is used in the following for creating the time-dependent solution.

You may have noticed the four —- after some of the comments. These are used for creating headlines in the script making it easier to navigate. Use the button “show document outline” in the top right corner of the editor pane (five vertical lines, three are indented).

The following lines are used to create labels for the ranges to be used later in the process. This in order to reflect the changes in parameters and ranges without rewriting.

#plotting out....-----
out_df <- as.data.frame(out1)
#wide to long ----
out_df_long <- out_df %>%
    pivot_longer(
        c(2:5),
        names_to = "Hormones",
        values_to = "Value",
        values_drop_na = TRUE
    )
#create labels for the ranges----
parRanges <- data.frame(min = c(s1min=s1min, 
                                k134min=k134min, 
                                k21min=k21min,
                                c2min=c2min,
                                k234min=k234min,
                                k32min=k32min, 
                                k43min=k43min,
                                k42min=k42min,
                                k423min=k423min 
                                
),
max = c(s1max=s1max, 
         k134max=k134max, 
         k21max=k21max,
         c2max=c2max,
         k234max=k234max,
         k32max=k32max, 
         k43max=k43max,
         k42max=k42max,
         k423max=k423max 
));
#rownames for labels
l_rownames <-c("s1",
               "k134",
               "k21",
               "c2",
               "k234",
               "k32",
               "k43", 
               "k42",
               "k423"
)

#before you attach rownames!!
#To create labels for the plots dynamically
labels <-data.frame(parRanges)
labels$parm <-l_rownames

Highlight the last points

Highlight the last points of the time dependent solution, by selecting the last point.

#prepare for last points hormones and value //TIME
Last10_3 <- out_df %>%
    filter(time == 20)
#Pivot_longer end points //TIME----
Last10_3 <-Last10_3 %>%
    pivot_longer(
        c(2:5),
        names_to = "Hormones",
        values_to = "Value",
        values_drop_na = TRUE
    )

The equilibrium solution

This chunk is part of the script from the supplemental material.

# function to find the equilibrium solution and vary selected parameters 

fs2= function(p) { 
    s1 <- steady(y=y0,func=fn1,parms=p,method="runsteady") 
    return(data.frame(TRH=s1$y[1],TSH=s1$y[2],FT4=s1$y[3],FT3=s1$y[4])) }

# vary k32 within a range 
parRanges <- data.frame(min = c(s1min=s1min, 
                                k134min=k134min, 
                                k21min=k21min,
                                c2min=c2min,
                                k234min=k234min,
                                k32min=k32min, 
                                k43min=k43min,
                                k42min=k42min,
                                k423min=k423min
                                ),  
                        max = c(s1max=s1max, 
                                k134max=k134max, 
                                k21max=k21max,
                                c2max=c2max,
                                k234max=k234max,
                                k32max=k32max, 
                                k43max=k43max,
                                k42max=k42max,
                                k423max=k423max
                                )); 
rownames(parRanges) <- c("s1",
                         "k134",
                         "k21",
                         "c2",
                         "k234",
                         "k32",
                         "k43", 
                         "k42",
                         "k423"); 

# find equilibrium solutions 
r1 <- sensRange(func=fs2,parms = k,parRange=parRanges, 
                sensvar=c("FT4","FT3","TSH","TRH"), num=1000);

Change r1 to dataframe

Reason: To be able to use the values in the plotting process

#plotting out....-----
r1_df <- as.data.frame(r1)

Selected points

Here we select points in r1_df very close to one. This can be changed if you want to highligt other values. The value (near one) changes a little - be aware (if precision matters). In line two the lowest displayed value are selected.

#prepare for last points hormones and value //EQUIL
#rm(data_k32_min)
data_k32 = subset(r1_df, c(k32 > 0.999 & k32 < 1.1))
data_k32_min = data_k32 %>% slice_min(data_k32$k32)

data_k32_min <-data_k32_min %>%
    pivot_longer(
        c(10:13),
        names_to = "Hormones",
        values_to = "Value",
        values_drop_na = TRUE
    )

Collect the two plots in one

When using only range for k32 and not any other - the figures in the two plots should be close to equal.

eq_time=
    time_10_3 +
    eq_10_3 +
    plot_annotation(tag_levels = 'A') +
    plot_layout(ncol = 2, byrow=TRUE)
eq_time