R assignment Day 1 for BIOS 20151 (summer 2025)

Work through R tutorials 1 and 2 prior to beginning this assignment

You will learn these skills

• Work with arithmetic operators
• Assign values to variables
• Print out variable values
• Assign vector variables
• Make plots using vectors
• Overlay plots and add legend

Grading and submitting

Type your name in the author field at the top (please don’t change anything else in that section), edit the code chunks in this document and type answers after the questions that ask for verbal answers as text outside the code chunks. You will be graded on correct completion of each task in each part, including both coding and verbal questions. There are numeric answers provided so you can check that your code is correct.

IMPORTANT: Before submitting, click the Render button at the top to make sure all the code runs from top to bottom. Submit the resulting HTML report through the link on Canvas.

Part 1: Blood volume pumped by heart

Take the model of a heart pumping with a constant rate, where V is the total volume of blood pumped by the heart over time, S is the stroke volume, R is the heart rate and t is time:

\[ V = SRt \]

1.1 For a heart beating at 80 beats per minute with stroke volume of 80 mL over 80 minutes, assign variables for the stroke volume, heart rate, and time (come up with your own descriptive names). Calculate the total volume pumped by the heart by turning the equation into a line of R code, assign the result to a variable called Vol and print out the result. ANSWER: 512000 (mL)

stroke_volume <- 80
heart_rate <- 80
time_minutes <- 80

Vol <- stroke_volume * heart_rate * time_minutes

print(Vol)

[1] 512000

1.2 Copy your script from above and modify the line that assigns the time variable to instead assign it a vector from 0 to 1440 minutes (24 hours) and keep the same values for stroke volume and heart rate. Then your calculation for total volume should produce a vector of the same length as time. To check your calculation, print out the value of the vector after 1440 minutes (corresponds to index 1441, since index 1 corresponds to 0 minutes). ANSWER: 9216000 (mL).

stroke_volume <- 80
heart_rate <- 80
time_minutes <- seq(0, 1440)

Vol <- stroke_volume * heart_rate * time_minutes[1441]

print(Vol)

[1] 9216000

1.3 Change the heart rate to 50 beats per minute and re-calculate the total blood volume with the new parameter, assigning it to a new variable (e.g. Vol2). How long does it take for the heart at this heart rate to pump 1000 liters (1 million mL) of blood? Try different indices and print out the one that gives the total volume closest to the target.

stroke_volume <- 80
heart_rate <- 50

time_minutes <- 250

Vol2 <- stroke_volume * heart_rate * time_minutes

print(time_minutes)

[1] 250

1.4 Plot the first volume variable (with heart rate of 80) vs time using the plot function with type line with black line color and label your axes. Plot the second volume variable (with heart rate of 50) on top of the first one using the lines() function with a different line color and add a legend that specifies the two different parameter values with their line colors.

stroke_volume <- 80
heart_rate <- 80

time <- seq(0, 80)
volume = stroke_volume * heart_rate * time

heart_rate <- 80
plot(time, volume, xlab = "Time (minutes)", ylab = "Volume (mL)", main = "Volume Blood Pumped over Time", type = 'l')

heart_rate_2 <- 50
volume_2 = stroke_volume * heart_rate_2  * time
lines(time, volume_2 , col = 'red')

legend("bottomright", legend = c('Heart Rate = 80bpm', 'Heart Rate = 50bpm'), col = c('black', 'red'), lty=1)

1.5 Describe the type of function and report the effect that the heart rate parameter has on its graph.

This is a linear function, and the heart rate parameter changes its slope.

Part 2: Pharmacokinetics

Many drugs are metabolized (transformed into easily removed compounds) in the bloodstream with a rate proportional to their concentration. This can be modeled using an exponential function with rate constant \(k\), starting with initial concentration \(C_0\):

\[ C = C_0 e^{-kt} \]

2.1 Start with a peak concentration of \(C_0=40\) (mcg/mL) and the metabolic decay rate of \(k=0.27\) (1/hour) (the decay rate for acetaminophen, also known as Tylenol) calculate the concentration of the drug after 12 hours, by assigning values to the variables on the right side of the equation, performing the calculation and assigning the result to a variable; then print out the result. ANSWER: about 1.6 (mcg/mL)

C_0 <- 40
k <- 0.27

t <- 12

acet_conc_12 = C_0 * exp(1) ^ (-k * t)


print(acet_conc_12)

[1] 1.566556

2.2 Copy your code from the chunk above and modify the line that assigns the time variable to instead assign it a vector from 0 to 12 hours; keep the same values for all the other variables. Then your calculation for acetaminophen concentration should produce a vector of the same length as time. To check the calculation, print out the value of concentration after 12 hours (index 13 of the concentration vector) ANSWER: about 1.6 (mcg/mL). How long does it take for the concentration to be halved (reach 20 mcg/mL)? Try different indices and print out the one that gives the total volume closest to the target.

C_0 <- 40
k <- 0.27

t <- seq(0, 12)

acet_conc_12 = C_0 * exp(1) ^ (-k * t)


print(acet_conc_12[13])

[1] 1.566556

print(3)

[1] 3

2.3 Change the metabolic decay rate to 0.17 (the value for rifampicin, a drug used to treat tuberculosis), keeping the other values the same, recalculate the concentration vector with the new parameter, and assign it to a new variable. To check the calculation, print out the value of concentration after 12 hours (index 13 of the concentration vector) ANSWER: about 5.2 (mcg/mL). How long does it take for the concentration to be halved (reach 20 mcg/mL)? Try different indices and print out the one that gives the total volume closest to the target.

C_0 <- 40
k <- 0.17

t <- seq(0, 12)

rifa_conc_12 = C_0 * exp(1) ^ (-k * t)

print(rifa_conc_12[13])

[1] 5.201148

print(5)

[1] 5

2.4 Plot the concentration of acetaminophen vs time using the type line with black line color and label your axes. Overlay the plot of concentration of rifampicin on top of the previous one using the lines() function with a different color and add a legend that specifies the two different parameter values with their line colors..

t = seq(0, 12)

k = 0.27
C = C_0 * exp(1) ^ (-k * t)

plot(t, C, xlab = "Time (hours)", ylab = "Concentration (mcg/mL)", main= "Drug Concentration over Time", type = 'l')

k_2 = 0.17
C_2 = C_0 * exp(1) ^ (-k_2 * t)
lines(t, C_2, col = "red")

legend("topright", legend = c('Acetaminophen', 'Rifampicin'), col = c('black', 'red'), lty=1)

2.5 Based on the plots, describe the type of function and report the effect the decay rate parameter has on its graph. Report the differences in half-life you found for for acetominophen and rifampicin.

This function is exponential decay, and the effect of a lower drug concentration decay rate has on this graph is slowing down the rate of drug metabolism. Rifampicin decays at roughly half the speed of Acetominphen.

Part 3: Population dynamics

The logistic model for a population P that is limited to a certain number, called carrying capacity K (population size), can be written as follows, with t time (years), r the intrinsic growth rate (per year) and A a dimensionless constant.

\[ P= \frac{K}{1 + A e^{-rt}} \]

3.1 Let \(K=10000\), \(r=0.1\), and \(A = 40\) and assign these values to informatively named variables. Calculate the population size after 25 years by turning the equation above into a line of R code, assigning the result to a variable, and print out the result. ANSWER: about 2335 individuals

carr_capac = 10000
growth_rate = 0.1
A = 40

time_yrs = 25

pop_25_yrs = carr_capac / (1 + A * exp(1) ^ (- growth_rate * time_yrs))

print(pop_25_yrs)

[1] 2334.594

3.2 Copy your code from the chunk above and modify the line that assigns the time variable to instead assign it a vector from 0 to 100 years (one for each year); keep the same values for all the other variables. Then your calculation for the population should produce a vector of the same length as time. Check your calculation by printing out the value of the vector after 15 years (index 16 of the vector) ANSWER: about 1008 individuals. Print out the initial value of the population with these parameters (index 1 of the population vector).

carr_capac = 10000
growth_rate = 0.1
A = 40

time_yrs_range = seq(0, 100)

pop_15_yrs = carr_capac / (1 + A * exp(1) ^ (- growth_rate * time_yrs_range[16]))
print(pop_15_yrs)

[1] 1007.536

pop_1_yrs = carr_capac / (1 + A * exp(1) ^ (- growth_rate * time_yrs_range[1]))
print(pop_1_yrs)

[1] 243.9024

3.3 Change the parameter A to 4, keeping the other values the same, recalculate the population vector with the new parameter, and assign it to a new variable. Check your calculation by printing out the value of the vector after 15 years (index 16 of the vector) ANSWER: about 5284 individuals. Print out the initial value of the population with the new parameter.

carr_capac = 10000
growth_rate = 0.1
A = 4

time_yrs_range = seq(0, 100)

pop_15_yrs_capac_4 = carr_capac / (1 + A * exp(1) ^ (- growth_rate * time_yrs_range[16]))
print(pop_15_yrs_capac_4)

[1] 5283.958

pop_1_yrs_capac_4 = carr_capac / (1 + A * exp(1) ^ (- growth_rate * time_yrs_range[1]))
print(pop_1_yrs_capac_4)

[1] 2000

3.4 Plot the population (with A=40) vs time using the type line with black line color and label your axes. Overlay the plot of the population (with A=4) on top of the first one using the lines() function with red line color and add a legend that specifies the two different parameter values with their line colors.

time_yrs_range = seq(0, 100)

A = 40
pop = carr_capac / (1 + A * exp(1) ^ (- growth_rate * time_yrs_range))

plot(time_yrs_range, pop, xlab = "Time (years)", ylab = "Population (individuals)", main= "Population over Time", type = 'l')

A_2 = 4
pop_2 = carr_capac / (1 + A_2 * exp(1) ^ (- growth_rate * time_yrs_range))
lines(time_yrs_range, pop_2, col = "red")

legend("bottomright", legend = c('A = 40', 'A = 4'), col = c('black', 'red'), lty=1)

3.5 Based on the plots, describe the type of function and report the effect that parameter A has on its graph.

This is a logistic growth function, and a lower A had the effect of a higher population being able to be sustained earlier in the population’s existence.