Merge pull request #16 from cropsinsilico/SoybeanParameterization

Soybean parameterization

Author: yufenghe22

CMakeLists.txt

@@ -117,8 +117,7 @@ include_directories(${SUNDIALS_INCLUDE_DIRS})
 include_directories("${PROJECT_SOURCE_DIR}/include")
 
 target_link_libraries(EPhotosynthesis ${SUNDIALS_LIBRARIES})
-
-target_link_libraries(EPhotosynthesis boost_regex)
+target_link_libraries(EPhotosynthesis ${Boost_LIBRARIES})
 target_link_libraries(ePhoto EPhotosynthesis)
 
 install(TARGETS EPhotosynthesis DESTINATION lib)

README.md

@@ -1,5 +1,9 @@
 # ePhotosynthesis
 
+On biocluster (Linux), I manage to build it successfully using a conda environment. First, initialize a new conda env so that the package path is on your local drive with r/w permissons. Second, install both boost and sundial through conda. Third, follow the following steps for cmake and make. This way, the cmake should be able to find the boost and sundial automatically. Of course, be sure to activate the related env the next time after login.
+
+Note: The sundials version I used is 5.7.0. The newer one 6.6.1 is not compatible!
+
 This is a C++ port of the ePhotosynthesis Matlab model code. The code is built with
 *CMake*.
 
@@ -42,4 +46,4 @@ The ePhotosynthesis executable takes the following arguments:
  -b,--begintime The starting time for the calculations (default is 0)
  -s,--stoptime  The ending time for the calculations (default is 5000)
  -z,--stepsize  The step size for the calculations (default is 1)
-```
+```

bin/ePhotosynthesis.cpp

@@ -126,7 +126,7 @@ int main(int argc, const char* argv[]) {
         theVars->TestATPCost = stoi(inputs.at("ATPCost"), nullptr);
         theVars->record = record;
         theVars->useC3 = useC3;
-        theVars->RUBISCOMETHOD = 1;
+        theVars->RUBISCOMETHOD = 2;
         modules::PR::setRUBISCOTOTAL(3);
         if (debugDelta)
             dbglvl += 8;

include/Variables.hpp

@@ -106,9 +106,11 @@ class Variables {
     double Tp = 0;
     double alfa = 0.;
     double fc = 0.;
+    double Phi_max = 0.;
     double lightParam = 0.;
-    const double alpha1 = 1.205;
-    const double alpha2 = 2.06;
+    const double alpha1 = 0.87;
+    const double alpha2 = 1.03;
+    double sensitivity_sf = 1.0; //a scaling factor for enzymes
 
     // Parameters
     arr PR_Param = zeros(2);

src/Condition.cpp

@@ -28,6 +28,7 @@
 #include "globals.hpp"
 #include <boost/algorithm/string_regex.hpp>
 #include <boost/regex.hpp>
+#include <sstream>
 
 const boost::regex token("\\s+");
 

src/drivers/EPS_Drive.cpp

@@ -24,11 +24,16 @@
  *
  **********************************************************************************************************************************************/
 
+#include <math.h>
 #include "Variables.hpp"
 #include "globals.hpp"
 #include "drivers/EPS_Driver.hpp"
 #include "modules/EPS.hpp"
 #include "modules/PS.hpp"
+//#include <iostream>
+//#include <fstream>
+
+//using namespace std;
 
 using namespace ePhotosynthesis;
 using namespace ePhotosynthesis::modules;
@@ -49,10 +54,13 @@ void EPSDriver::setup() {
     SYSInitial(inputVars);
     //time = tglobal;
     inputVars->Tp = this->Tp;
-    inputVars->alfa = 0.85;
+    inputVars->Phi_max = 0.63;
     PS::setJmax(inputVars->EnzymeAct.at("Jmax"));
-    inputVars->fc = 0.15;
-    PS::setTheta(0.7);
+    inputVars->fc = 0.2;
+
+    double Tp = inputVars->Tp;
+    double Theta = 0.76 + 0.01713 * Tp - 3.75 * pow(Tp,2.0) / 10000.0;//Yufeng: match Farquhar Matlab
+    PS::setTheta(Theta);
     PS::setbeta(0.7519);
     inputVars->BF_FI_com = true;
     inputVars->EnzymeAct.at("V1") *= inputVars->alpha1;
@@ -67,6 +75,11 @@ void EPSDriver::setup() {
     inputVars->EnzymeAct.at("V13") *= inputVars->alpha2;
     inputVars->EnzymeAct.at("V23") *= inputVars->alpha2;
 
+    //we scale up some enzymes 
+    inputVars->EnzymeAct.at("V2") *= inputVars->sensitivity_sf;
+    inputVars->EnzymeAct.at("V13") *= inputVars->sensitivity_sf;
+    //
+
     inputVars->PR_PS_com = true;     // This is a variable indicating whether the PR model is actually need to be combined with PS or not. If 1 then means combined; 0 means not.
 
 
@@ -110,10 +123,17 @@ void EPSDriver::setup() {
 void EPSDriver::getResults() {
     EPSCondition* eps_int_con = new EPSCondition(intermediateRes);
     arr temp = EPS::MB(time, eps_int_con, inputVars);
-    results = zeros(1);
+    results = zeros(3);
     const double Arate = (inputVars->PS_Vel.v1 - inputVars->PR_Vel.v131) * inputVars->AVR;
     delete eps_int_con;
+
+//    ofstream myfile;
+//    myfile.open ("hello.txt");
+//    myfile << "Writing this to a file.\n";
+//    myfile.close();
     results[0] = Arate;
+    results[1] = inputVars->PS_Vel.v1 * inputVars->AVR; //carboxylation
+    results[2] = inputVars->PR_Vel.v131 * inputVars->AVR; //photorespiration 
 }
 
 EPSCondition* EPSDriver::EPS_Init() {

src/ini/BF_Ini.cpp

@@ -65,7 +65,7 @@ BFCondition* BF::_init(Variables *theVars) {
         // CPSi=1.0131;% 1.0237 WY201803
         // cNADPHsyn=1.094468408;%1.0388 WY201803
         if (theVars->lightParam == 0.) {
-            const double light_scaler = theVars->alfa * (1. - theVars->fc);
+            const double light_scaler = theVars->Phi_max * (1. - theVars->fc);
             theVars->lightParam = theVars->TestLi * 30. * light_scaler;
         }
 

src/ini/FI_Ini.cpp

@@ -53,8 +53,8 @@ FICondition* FI::_init(Variables *theVars) {
     if (theVars->useC3) {
         FI::cpsii = 1.;
         if (theVars->lightParam == 0.) {
-            const double light_scaler = theVars->alfa * (1 - theVars->fc);
-            theVars->lightParam = theVars->TestLi * 30 * light_scaler;
+            const double light_scaler = theVars->Phi_max * (1. - theVars->fc);
+            theVars->lightParam = theVars->TestLi * 30. * light_scaler;
         }
         theVars->FI_RC.kA_d = theVars->EnzymeAct.at("kA_d");         // The rate constant of heat dissipation from peripheral antenna Lazar (1999), 0.25~1 *10^(9)
         theVars->FI_RC.kA_f = theVars->EnzymeAct.at("kA_f");         // The rate constant of fluorescence emission from peripheral antenna Lazar 1999, with a lifetime of 5 ns at closed reaction center

src/ini/PS_Ini.cpp

@@ -27,11 +27,12 @@
 #include "Variables.hpp"
 #include "modules/PS.hpp"
 
-#define R 8.314
-#define c_c 38.28
-#define dHa_c 80.99
-#define c_o 14.68
-#define dHa_o 23.72
+//Yufeng: I changed these constants to match BioCro's
+#define R 8.314E-3 //KJ mole-1 K-1
+#define c_c 38.05
+#define dHa_c 79.43 
+#define c_o 20.30 
+#define dHa_o 36.38 
 #define RegFactor 1.
 
 using namespace ePhotosynthesis;
@@ -223,7 +224,7 @@ PSCondition* PS::_init(Variables *theVars) {
     PS_con->PenP = 0.25;
 
     if (theVars->useC3) {
-        PS::PS_C_CP = 22.;    //   Global constant for the total phosphate
+        PS::PS_C_CP = 30.;    //   Global constant for the total phosphate
         PS::PS_C_CA = 1.5;   //   Global constant for the total adenylates
         PS::PS_C_CN = 1.;    //   Global constant for the cytosolic Phosphate concentration;
         if (theVars->GRNC == 1 && theVars->CO2_cond > 0) {
@@ -248,8 +249,10 @@ PSCondition* PS::_init(Variables *theVars) {
 
         //PsKM11_0 = ;
         //PsKM12_0 = ;  // O2 1 RuBP+CO2->2PGA  0.28 DEFAUL.
-        PS::KM11 = 0.0097 * exp(c_c - dHa_c * 1000. / (R * (theVars->Tp + 273.15))) / 272.38;
-        PS::KM12 = 0.244 * exp(c_o - dHa_o * 1000. / (R * (theVars->Tp + 273.15))) / 165.82;
+        double kc_25 = exp(c_c - dHa_c / (R * (25.0 + 273.15)));
+        double ko_25 = exp(c_o - dHa_o / (R * (25.0 + 273.15))); 
+        PS::KM11 = 0.0097 * exp(c_c - dHa_c / (R * (theVars->Tp + 273.15))) / kc_25;
+        PS::KM12 = 0.244  * exp(c_o - dHa_o / (R * (theVars->Tp + 273.15))) / ko_25; 
 
         PS::KM13 = 0.02;   //  RuBP 1 RuBP+CO2->2PGA
         PS::KI11 = 0.84;   // PGA
@@ -373,7 +376,7 @@ PSCondition* PS::_init(Variables *theVars) {
         PS::PsV10 = PS::PsV10_0 * pow(Q10_10, (theVars->Tp - 25.) / 10.);
         PS::PsV13= PS::PsV13_0 * pow(Q10_13, (theVars->Tp - 25.) / 10.);
         PS::PsV23 = PS::PsV23_0 * pow(Q10_23, (theVars->Tp - 25.) / 10.);
-        PS::I2 = theVars->TestLi * theVars->alfa * (1. - theVars->fc) / 2.;
+        PS::I2 = theVars->TestLi * theVars->Phi_max * (1. - theVars->fc) / 2.;
         PS::J = (I2 + PS::Jmax - sqrt(pow(I2 + PS::Jmax, 2.) - 4. * PS::Theta * I2 * PS::Jmax)) / (2. * PS::Theta);
 
     } else {

src/rate/PR_Rate.cpp

@@ -80,8 +80,11 @@ void PR::_Rate(const double t, const PRCondition* const PR_con, Variables *theVa
             theVars->PR_Vel.v111 = PrV111t * theVars->O2_cond / (theVars->O2_cond + PR::KO *
                                                                  (1. + theVars->CO2_cond / PR::KC));
 
-            if (RuBP < PR::RUBISCOTOTAL)
-                theVars->PR_Vel.v111 = theVars->PR_Vel.v111 * RuBP / PR::RUBISCOTOTAL;
+//            if (RuBP < PR::RUBISCOTOTAL)
+//                theVars->PR_Vel.v111 = theVars->PR_Vel.v111 * RuBP / PR::RUBISCOTOTAL;
+            if (RuBP < PS::getPsV1() / 2.) {
+                theVars->PR_Vel.v111 = theVars->PR_Vel.v111 * RuBP / (PS::getPsV1()/ 2.);
+            }
 
         } else if (theVars->RUBISCOMETHOD == 1) {
             theVars->PR_Vel.v111 = PrV111 * theVars->O2_cond / (theVars->O2_cond + PR::KO *