#Part 3 
K<- .1 #torr
pO2 <- seq(0,105,.01) #mmHg
Y <-((pO2*K)+(3*(pO2^2)*(K^2))+(3*(pO2^3)*(K^3))+((pO2^4)*(K^4)))/(1+(4*pO2*K)+(6*(pO2^2)*(K^2))+(4*(pO2^3)*(K^3))+((pO2^4)*(K^4)))*100
#assign x range input and y output

ExperimentalData = read.csv("C:/Users/maria/Documents/Physio/Coding Labs/Experimental_Data.csv") 
#fetch data frame


plot(pO2,Y, xlab= "pO2 (mmHg)", ylab= "Hemoglobin Saturation (%)", main= "Hemoglobin Oxygen Dissociation Curve (Independent Binding Simulation + Experimental Data)", xlim= c(0,100), ylim= c(0,100), type="l", col= c(1)) 
lines(ExperimentalData, type= "l", col= c(2))
#plot independent binding and add experimental data curve

legend("bottomright",legend = c("independent binding","experimental data"),col=c(1,2), lwd=1, cex=1)

#add legend

#1)the independent binding curve is hyperbolic, shows greater percentage of hemoglobin saturated at a lower pO2 value, lower p50. the experimental binding curve is sigmoidal, shows cooperative binding occurs, saturation point for O2 occurs at higher pO2, higher p50. 

#2)

plot(0,0, xlab= "pO2 (mmHg)", xlim= c(0,100), ylim= c(0,100), ylab= "Hemoglobin Saturation (%)", main= "Hemoglobin Oxygen Dissociation Curve (Independent Binding Simulation w/ 3 K test values + Experimental Data)") #base plot
lines(ExperimentalData, type= "l", col=c(4)) #experimental data curve
abline(h=50, col= c(5)) #horizontal line @ 50

KRange <- c(.05, .025, .0125) #torr
for(i in 1:4){
  YRange <-((pO2*KRange[i])+(3*(pO2^2)*(KRange[i]^2))+(3*(pO2^3)*(KRange[i]^3))+((pO2^4)*(KRange[i]^4)))/(1+(4*pO2*KRange[i])+(6*(pO2^2)*(KRange[i]^2))+(4*(pO2^3)*(KRange[i]^3))+((pO2^4)*(KRange[i]^4)))*100
  lines(pO2,YRange, col=c(i))
}
#3 K values 

legend("topleft", legend= c(KRange, "Experimental Data"), fill= c(1:4))

#3) I guess K at 1/40 is the best, but independent binding does not represent hemoglobin binding. None of the three K values is sigmoidal. Even the best plot is hyperbolic.
#Part 4

pO2 <- seq(0,105,.01) #mmHg
p50 <- 47.5 #mmHg
Y <- ((pO2^4)/((pO2^4)+ (p50^4)))*100 #%

plot(pO2,Y, xlab= "pO2 (mmHg)", ylab= "Hemoglobin Saturation (%)", main= "Hemoglobin Oxygen Dissociation Curve (Simultaneous Binding Simulation Data + Experimental Data)", xlim= c(0,100), ylim= c(0,100), type="l", col= c(1)) #plot simultaneous binding
lines(ExperimentalData, type= "l", col= c(2)) #add Experimental Data

legend("bottomright",legend = c("simultaneous binding","experimental data"),col=c(1,2), lwd=1, cex=1)

#4) the two curves are pretty much the same, I toggled with the p50 a little, but 40-50 mmHg gives a curve slightly less to slightly more sigmoidal than the experimental curve

#5)
plot(0,0, xlab= "pO2 (mmHg)", ylab= "Hemoglobin Saturation (%)", main= "Hemoglobin Oxygen Dissociation Curve (Simultaneous Binding Simulation Data + Experimental Data)", xlim= c(0,100), ylim= c(0,100), type="l") #baseplot
lines(ExperimentalData, col= c(6)) #add Experimental Data Curve
abline(h=50)

p50Range <- c(43.5 ,45.5,  47.5, 49.5, 51.5)
for (i in 1:5){
  YRange <- ((pO2^4)/((pO2^4)+ (p50Range[i]^4)))*100
  lines(pO2,YRange, col= c(i))
} #add p50 test values

legend("topleft", legend= c(p50Range, "Experimental Data"), fill= c(1:6), lwd= 1, cex= 1 )

#6) best fit is 47.5. It's sigmoidal. Because the plort matches the sigmoidal so well, I think O2 binding in hemoglobin could be simultaneous, but none of the curves matched exactly.
#Part 5

pO2 <- seq(0,105,.01) #mmHg
p50 <- 47.5 #mmHg
n <- 4
Y <- ((pO2^n)/((pO2^n)+ (p50^n)))*100 #%

plot(pO2,Y, xlab= "pO2 (mmHg)", ylab= "Hemoglobin Saturation (%)", main= "Hemoglobin Oxygen Dissociation Curve (Cooperative Binding Simulation Data + Experimental Data)", xlim= c(0,100), ylim= c(0,100), type="l", col= c(1)) #plot cooperative binding
lines(ExperimentalData, type= "l", col= c(2)) #add Experimental Data

legend("bottomright",legend = c("cooperative binding","experimental data"),col=c(1,2), lwd=1, cex=1)

#7) it's the same as with simultaneous binding. The two curves are both sigmoidal, but the simulation curve increases with a steeper slope

plot(0,0, xlab= "pO2 (mmHg)", ylab= "Hemoglobin Saturation (%)", main= "Hemoglobin Oxygen Dissociation Curve (Cooperative Binding Simulation Data + Experimental Data)", xlim= c(0,100), ylim= c(0,100), type="l") #baseplot
lines(ExperimentalData, col= c(5)) #add Experimental Data Curve
abline(h=50) # horizontal line

n <- c(3.3 ,3.7, 4.1, 4.5)
for (i in 1:4){
  YRange <- ((pO2^n[i])/((pO2^n[i])+ (p50^n[i])))*100
  lines(pO2,YRange, col= c(i))
} #add n test values

legend("topleft", legend= c(n, "Experimental Data"), fill= c(1:5), lwd= 1, cex= 1 )

#make legend

#8) best n seems to be 3.7

#9) shape of all plots is sigmoidal. I think O2 binding is cooperative because you can adjust the simultaneous plot even more to get it closer to the experimental data curve
#Part 6

#increasing n makes steeper sigmoidal slope
#increasing p50 shifts right

#10)

tdata <-read.csv("C:/Users/maria/Documents/Physio/Coding Labs/Temp_Data_Exp.csv") 
#fetch data frame

pO2 <- seq(0,105,.01) #mmHg

#32 degrees
n_32 <- 3.4
p50_32 <- 26
Y_32 <- ((pO2^n_32)/((pO2^n_32)+(p50_32^n_32)))*100 #%
plot( tdata [,1], tdata [,2], xlab= "pO2 mmHg", xlim= c(0,100), ylim= c(0,100), ylab= "% Saturation", main= " Hemoglobin Dissociation Curve 32 degrees Celsius", type= "l")
lines(pO2, Y_32, type="l", col=c(2))
legend("topleft", legend= c("Experimental Data", "Simulation Data"), col= c(1,2), lwd= 1, cex= 1)

#35 degrees
n_35 <- 3.4
p50_35 <- 39
Y_35 <- ((pO2^n_35)/((pO2^n_35)+(p50_35^n_35)))*100 #%
plot( tdata [,3], tdata [,4], xlab= "pO2 mmHg", xlim= c(0,100), ylim= c(0,100),ylab= "% Saturation", main= " Hemoglobin Dissociation Curve 35 degrees Celsius", type= "l")
lines(pO2, Y_35, type="l", col=c(2))
legend("topleft", legend= c("Experimental Data", "Simulation Data"), col= c(1,2), lwd= 1, cex= 1)

#38 degrees
n_38<- 3.4
p50_38 <- 42
Y_38 <- ((pO2^n_38)/((pO2^n_38)+(p50_38^n_38)))*100 #%
plot( tdata [,5], tdata [,6], xlab= "pO2 mmHg", xlim= c(0,100), ylim= c(0,100),ylab= "% Saturation", main= " Hemoglobin Dissociation Curve 38 degrees Celsius", type= "l")
lines(pO2, Y_38, type="l", col=c(2))
legend("topleft", legend= c("Experimental Data", "Simulation Data"), col= c(1,2), lwd= 1, cex= 1)

#41 degrees
n_41 <- 3.3
p50_41 <- 46
Y_41 <- ((pO2^n_41)/((pO2^n_41)+(p50_41^n_41)))*100 #%
plot( tdata [,7], tdata [,8], xlab= "pO2 mmHg", xlim= c(0,100), ylim= c(0,100),ylab= "% Saturation", main= " Hemoglobin Dissociation Curve 41 degrees Celsius", type= "l")
lines(pO2, Y_41, type="l", col=c(2))
legend("topleft", legend= c("Experimental Data", "Simulation Data"), col= c(1,2), lwd= 1, cex= 1)

#44 degrees
n_44 <- 3.4
p50_44 <- 49
Y_44 <- ((pO2^n_44)/((pO2^n_44)+(p50_44^n_44)))*100 #%
plot( tdata [,9], tdata [,10], xlab= "pO2 mmHg", xlim= c(0,100), ylim= c(0,100),ylab= "% Saturation", main= " Hemoglobin Dissociation Curve 44 degrees Celsius", type= "l")
lines(pO2, Y_44, type="l", col=c(2))
legend("topleft", legend= c("Experimental Data", "Simulation Data"), col= c(1,2), lwd= 1, cex= 1)

#47 degrees
n_47 <- 5
p50_47 <- 49
Y_47 <- ((pO2^n_47)/((pO2^n_47)+(p50_47^n_47)))*100 #%
plot( tdata [,11], tdata [,12], xlab= "pO2 mmHg",xlim= c(0,100), ylim= c(0,100), ylab= "% Saturation", main= " Hemoglobin Dissociation Curve 47 degrees Celsius", type= "l")
lines(pO2, Y_47, type="l", col=c(2))
legend("topleft", legend= c("Experimental Data", "Simulation Data"), col= c(1,2), lwd= 1, cex= 1)

#52 degrees
n_52 <- 2.15
p50_52 <- 54
Y_52 <- ((pO2^n_52)/((pO2^n_52)+(p50_52^n_52)))*100 #%
plot( tdata [,13], tdata [,14], xlab= "pO2 mmHg",xlim= c(0,100), ylim= c(0,100), ylab= "% Saturation", main= " Hemoglobin Dissociation Curve 52 degrees Celsius", type= "l")
lines(pO2, Y_52, type="l", col=c(2))
legend("topleft", legend= c("Experimental Data", "Simulation Data"), col= c(1,2), lwd= 1, cex= 1)

#11)
nvector <- c(3.4,3.4,3.4,3.3,3.4,5,2.15)
p50vector<- c(26,39,42,46,49,49,54)
tempvector <- c(32,35,38,41,44,47,52)

#12)
plot(tempvector,p50vector, xlab= "temperature (degrees Celsius)", ylab= "p50 (mmHg)", main= "p50 (mmHg) vs Temperature (Celsius)", type= "o")

#p50 vs T

plot(tempvector,nvector, xlab= "temperature (degrees Celsius)", ylab= "n", main= "n vs Temperature (Celsius)", type= "o")

#n vs T

#13)
p50linearfit <- lm(p50vector ~ tempvector)
print(p50linearfit)
## 
## Call:
## lm(formula = p50vector ~ tempvector)
## 
## Coefficients:
## (Intercept)   tempvector  
##      -6.558        1.214
cor(p50vector,tempvector)
## [1] 0.9217418
#p50= 1.214(temp) -6.558 mmHg; (slope 1.214, y intercept -6.558)
#correlation = .922

nlinearfit <- lm(nvector ~ tempvector)
print(nlinearfit)
## 
## Call:
## lm(formula = nvector ~ tempvector)
## 
## Coefficients:
## (Intercept)   tempvector  
##     4.03375     -0.01449
cor(nvector,tempvector)
## [1] -0.1217746
#n= 4.03375(temp) - .01449; (slope 4.03375, y int -.01449)
#correlation= -.122

#p50 vs temp is lineaer, as temp increases, so does p50 val
#n is not linear, looking at the graph it stays pretty constant, there are two outliers that skew the data at the higher temperatures. the sigmoidal nature of the hemoglobin oxygen dissociation curve shouldn't change

#14)
par(mar= c(4,4,1,4))
plot(tempvector, p50vector, xlab= "Temperature (degrees Celsius)", xlim= c(30,55), ylim= c(20,55), ylab= "p50 values (mmHg)", main= "p50 and n vs Temperature", type= "l", col= c(1))
abline(p50linearfit, col= c(3))
par(new=TRUE)
plot(tempvector, nvector, xaxt= "n", yaxt= "n", ylab= "", xlab= "", col= c(2), type= "l", lty= 2)
abline(nlinearfit,col= c(4))
axis(side= 4)
mtext("n values", side=4, line= 3)
legend("bottom", c("p50 v temp", "n v temp", "p50 v temp linear fit", "n v temp linear fit"), col= c(1,2, 3,4), lty= c(1,2), cex=.5)

#Part 7
pHvector <- c(7.61,7.36,7.15,6.92)
p50vector <- c(19.1,23.4,30.9,40.7)

#15)
plot(p50vector, pHvector, xlab= "p50 (mmHg)", ylab= "pH", main= "pH vs p50")

#16)
pHvp50fit <- lm(pHvector ~ p50vector)
pHvp50fit
## 
## Call:
## lm(formula = pHvector ~ p50vector)
## 
## Coefficients:
## (Intercept)    p50vector  
##     8.13091     -0.03053
cor(pHvector,p50vector)
## [1] -0.9816245
#slope= -.03053, y-int= 8.13091
#correlation= -.982

plot(p50vector,pHvector, xlab= "p50(mmHg)", ylab= "pH", main= "pH vs p50")
abline(pHvp50fit)

#it's an average fit

#17)

p50function <- function(pHvector){
  p50 = c(0,0,0,0)
  for (i in 1:4){
    p50[i] = ((pHvector[i]-8.13091)/(-.03053))
  }
  return(p50)
}
p50output <- p50function(pHvector)
p50output
## [1] 17.06223 25.25090 32.12938 39.66295
#17.06223, 25.25090, 32.12938,39.66295

#18)
pO2 <- seq(0,150,.01)
plot(0, 0, xlab= "pO2 (mmHg)", xlim= c(0,150), ylim= c(0,100), ylab= "Hemoglobin Saturation %", main= "Hemoglobin Dissociation Curve" )
p50 <- c(17.06223, 25.25090, 32.12938,39.66295)
for (i in 1:4){
  Hill <- ((pO2^2.8)/((pO2^2.8)+(p50[i]^2.8)))*100
  lines(pO2,Hill, col= c(i))
}
legend("bottomright", c("pH 7.61", "ph 7.36","pH 7.15","pH 6.92"), fill= c(1:4), lwd= 1, cex=1)

#19) As pH increases, the slope gets steeper. This means it reaches its p50 and becomes fully saturated at lower pO2 with higher pH.pH seems to increase affinity.

#20)
abline(v=40)
abline(v=100)

#21)
#pH has more ofan effect at 40 torr. This means that hemoglobin of arterial blood is not likely to unload with small pO2 changes no matter the pH. Venous blood is likely to unload with small changes in pO2, with higher pH making the changes more drastic.

#22)
Ymyo <- ((pO2)/(pO2 + 2.4))*100
lines(pO2, Ymyo, col=c(5))

#23)
#myoglobin binds more tightly (higher binding affinity) and is less likely to let go of its oxygen unless pO2 is really low, that way if there's no oxygen present, myoglobin can deliver the O2 to the muscles. Hemoglobin has a lower binding affinity, so it can release its oxygen more and distribute it when there are more normal levels of O2