SciCodePile/SciCode-Domain-Code
DATA1: Domain-Specific Code Dataset Dataset Overview DATA1 is a large-scale domain-specific code dataset focusing on code samples from interdisciplinary fields such as biology, chemistry, materials science, and related areas. The dataset is collected and organized from GitHub repositories, covering 178 different domain topics with over 1.1 billion lines of code. Dataset Statistics Total Datasets: 178 CSV files Total Data Size: ~115 GB Total Lines… See the full description on the dataset page: https://huggingface.co/datasets/SciCodePile/SciCode-Domain-Code.
42.4k
1"keyword","repo_name","file_path","file_extension","file_size","line_count","content","language"
2"Inhibition","simonetacannodelas/BioVL-Library","Rx_Fermentation_Ecoli_Glucose_Aerobic_Fedbatch.py",".py","11040","292","# -*- coding: utf-8 -*-3""""""4Created on Friday Jan 24 13:34:32 20205 6@author: simoca7""""""8from scipy.integrate import odeint 9#Package for plotting 10import math 11#Package for the use of vectors and matrix 12import numpy as np13import pandas as pd14import array as arr15from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas16from matplotlib.figure import Figure17import sys18import os19import matplotlib.pyplot as plt20from matplotlib.ticker import FormatStrFormatter21import glob22from random import sample23import random24import time25import plotly26import plotly.graph_objs as go27import json28from plotly.subplots import make_subplots29 30 31class Ecoli_Aero:32 def __init__(self, Control = False):33 self.Kap=0.5088 #g/L34 self.Ksa=0.012835 self.Kia=1.2602 #g/L36 self.Ks= 0.0381 #g/L37 self.Kis = 1.8383 #g/L affinity constant38 self.Ko = 0.000139 self.qAcmax= 0.1148 #g/gh40 self.qm=0.0133 #g/gh41 self.qOmax= 13.4*31.9988/1000 #g/gh42 self.qSmax= 0.635 #g/gh43 self.Yas= 0.8938 #g/g44 self.Yoa= 0.5221 #g/g45 self.Yos= 1.5722 #g/g46 self.Yxa= 0.5794 #g/g47 self.Yem = 0.5321 #g/g48 self.Yxsof = 0.229 # g/g49 self.pAmax= 0.2286 #gA/gXh50 self.kla= 22051 self.H= 1400 #Henry constant52 53 self.G0 = 4.9454 self.A0 = 0.012955 self.O0 = 9856 self.X0 = 0.1757 self.tau= 35 #response time58 self.t_start = 059 self.V0 = 260 self.F0 = 061 self.SFR = 1.562 self.t_expfb_start = 063 self.t_constfb_start = 1.564 self.t_end = 1565 66 #parameters for control, default every 1/24 hours:67 self.Control = Control68 self.coolingOn = True69 self.Contamination=False70 self.steps = (self.t_end - self.t_start)*2471 self.T0 = 3572 self.K_p = 2.31e+0173 self.K_i = 174 self.K_d = 075 self.Tset = 3076 self.u_max = 15077 self.u_min = 078 79 def update_param_value(self, param, new_value):80 81 class_name = str(self.__class__).split('.')[1].split(""'"")[0]82 param_current_value = self.__dict__.get(param, None)83 84 if param_current_value is None:85 print(""Class {} does not contain param with name {}"".format(class_name, param))86 return87 88 self.__dict__[param] = new_value89 print(""Value of param {} in class {} updated to {}"".format(param, class_name, new_value))90 91 92 def rxn(self, C,t , u, fc):93 #when there is no control, k has no effect94 k=195 #when cooling is off than u = 096 if self.coolingOn == False:97 u = 098 if self.Contamination == True:99 fc=np.random.randint(0,10)100 fc=fc/17101 102 if self.Control == True :103 #Cardinal temperature model with inflection: Salvado et al 2011 ""Temperature Adaptation Markedly Determines Evolution within the Genus Saccharomyces""104 #Strain E.coli W310105 Topt = 35106 Tmax = 45.48107 Tmin = 10108 T = C[5]109 if T < Tmin or T > Tmax:110 k = 0111 else:112 D = (T-Tmax)*(T-Tmin)**2113 E = (Topt-Tmin)*((Topt-Tmin)*(T-Topt)-(Topt-Tmax)*(Topt+Tmin-2*T))114 k = D/E115 # Volume balance116 if (t >= self.t_expfb_start):117 if (t < self.t_constfb_start):118 Fin = self.F0 * math.exp(self.SFR * (t - self.t_expfb_start))119 Fout = 0120 else:121 Fin = self.F0 * math.exp(self.SFR * (self.t_constfb_start - self.t_expfb_start))122 Fout = 0123 else:124 Fin = 0125 Fout = 0126 127 F = Fin - Fout128 129 qS = (self.qSmax/(1+C[1]/self.Kia))*(C[0]/(C[0]+self.Ks))130 qSof = self.pAmax*(qS/(qS+self.Kap))131 pA = qSof*self.Yas132 qSox = (qS-qSof)*(C[2]/(C[2]+self.Ko))133 qSan = (qSox-self.qm)*self.Yem*(0.488/0.391)134 qsA = (self.qAcmax/(1+(qS/self.Kis)))*(C[1]/(C[1]+self.Ksa))135 qA = pA - qsA136 mu = (qSox - self.qm)*self.Yem+ qsA*self.Yxa + qSof*self.Yxsof137 qO = self.Yos*(qSox-qSan)+qsA*self.Yoa138 139 #Solving the mass balances140 dGdt = (F/C[4])*(self.G0-C[0])-(qS*C[3])141 dAdt = qA*C[3]-((F/C[4])*C[1])142 dOdt = self.kla*(self.O0-C[2])-qO*C[3]*self.H143 dXdt = (mu - (F/C[4]))*C[3]144 dVdt = F145 if self.Control == True :146 '''147 dHrxn heat produced by cells estimated by yeast heat combustion coeficcient dhc0 = -21.2 kJ/g148 dHrxn = dGdt*V*dhc0(G)-dEdt*V*dhc0(E)-dXdt*V*dhc0(X)149 (when cooling is working) Q = - dHrxn -W ,150 dT = V[L] * 1000 g/L / 4.1868 [J/gK]*dE [kJ]*1000 J/KJ151 dhc0(EtOH) = -1366.8 kJ/gmol/46 g/gmol [KJ/g]152 dhc0(Glc) = -2805 kJ/gmol/180g/gmol [KJ/g]153 154 ''' 155 #Metabolic heat: [W]=[J/s], dhc0 from book ""Bioprocess Engineering Principles"" (Pauline M. Doran) : Appendix Table C.8 156 dHrxndt = dXdt*C[4]*(-21200) #[J/s] + dGdt*C[4]*(15580)- dEdt*C[4]*(29710) 157 #Shaft work 1 W/L1158 W = 1*C[4] #[J/S] negative because exothermic 159 #Cooling just an initial value (constant cooling to see what happens)160 #dQdt = -0.03*C[4]*(-21200) #[J/S] 161 #velocity of cooling water: u [m3/h] -->controlled by PID 162 163 #Mass flow cooling water164 M=u/3600*1000 #[kg/s]165 #Define Tin = 5 C, Tout=TReactor166 #heat capacity water = 4190 J/kgK167 Tin = 5168 #Estimate water at outlet same as Temp in reactor169 Tout = C[5]170 cpc = 4190171 #Calculate Q from Eq 9.47172 Q=-M*cpc*(Tout-Tin) # J/s 173 #Calculate Temperature change174 dTdt = -1*(dHrxndt - Q + W)/(C[4]*1000*4.1868) #[K/s]175 else: 176 dTdt = 0177 return [dGdt,dAdt,dOdt,dXdt,dVdt, dTdt]178 179 def solve(self):180 #solve normal:181 t = np.linspace(self.t_start, self.t_end, self.steps)182 if self.Control == False :183 u = 0184 fc= 1185 C0 = [self.G0, self.A0, self.O0, self.X0,self.V0,self.T0]186 C = odeint(self.rxn, C0, t, rtol = 1e-7, mxstep= 500000, args=(u,fc,))187 188 #solve for Control189 else:190 fc=0191 """"""192 PID Temperature Control:193 """"""194 # storage for recording values195 C = np.ones([len(t), 6]) 196 C0 = [self.G0, self.A0, self.O0, self.X0,self.V0,self.T0]197 self.ctrl_output = np.zeros(len(t)) # controller output198 e = np.zeros(len(t)) # error199 ie = np.zeros(len(t)) # integral of the error200 dpv = np.zeros(len(t)) # derivative of the pv201 P = np.zeros(len(t)) # proportional202 I = np.zeros(len(t)) # integral203 D = np.zeros(len(t)) # derivative204 205 for i in range(len(t)-1):206 #print(t[i])207 #PID control of cooling water208 dt = t[i+1]-t[i]209 #Error210 e[i] = C[i,5] - self.Tset 211 #print(e[i])212 if i >= 1:213 dpv[i] = (C[i,5]-C[i-1,5])/dt214 ie[i] = ie[i-1] + e[i]*dt215 P[i]=self.K_p*e[i]216 I[i]=self.K_i*ie[i]217 D[i]=self.K_d*dpv[i]218 219 self.ctrl_output[i]=P[i]+I[i]+D[i]220 u=self.ctrl_output[i]221 if u>self.u_max:222 u=self.u_max223 ie[i] = ie[i] - e[i]*dt # anti-reset windup224 if u < self.u_min:225 u =self.u_min226 ie[i] = ie[i] - e[i]*dt # anti-reset windup227 #time for solving ODE 228 ts = [t[i],t[i+1]]229 #disturbance230 #if self.t[i] > 5 and self.t[i] < 10:231 # u = 0 232 #solve ODE from last timepoint to new timepoint with old values 233 234 y = odeint(self.rxn, C0, ts, rtol = 1e-7, mxstep= 500000, args=(u,fc,))235 #update C0236 C0 = y[-1]237 #merge y to C238 C[i+1]=y[-1]239 return t, C240 def create_plot(self, t, C):241 figure = make_subplots(rows=2, cols=1)242 [self.G0, self.A0, self.O0, self.X0, self.V0, self.T0]243 G = C[:, 0]244 A = C[:, 1]245 B = C[:, 3]246 O = C[:, 2]247 V = C[:, 4]248 df = pd.DataFrame({'t': t, 'G': G, 'B': B, 'A': A, 'O': O, 'V': V})249 figure.append_trace(go.Scatter(x=df['t'], y=df['G'], name='Glucose'), row=1, col=1)250 figure.append_trace(go.Scatter(x=df['t'], y=df['O'], name='Oxygen'), row=1, col=1)251 figure.append_trace(go.Scatter(x=df['t'], y=df['B'], name='Biomass'), row=1, col=1)252 figure.append_trace(go.Scatter(x=df['t'], y=df['A'], name='Acetate'), row=1, col=1)253 # fig.update_layout(title=('Simulation of the model for the Scerevisiae'),254 # xaxis_title='time (h)',255 # yaxis_title='Concentration (g/L)')256 257 258 Tset = self.T0259 df2 = pd.DataFrame({'t': t, 'T': T, 'Tset': Tset})260 figure.append_trace(go.Scatter(x=df2['t'], y=df2['T'], name='Temperature'), row=2, col=1)261 figure.append_trace(go.Scatter(x=df2['t'], y=df2['Tset'], name='Set Value Temperature'), row=2, col=1)262 figure.update_layout(title=('Simulation of the model for the Saccharomyces cerevisiae'),263 xaxis_title='time (h)',264 yaxis_title='Concentration (g/L)')265 266 S = C[:, 0]267 A = C[:, 1]268 B = C[:, 3]269 df = pd.DataFrame({'t': t, 'Substrate': S, 'Biomass': B, 'Acetate': A})270 fig = go.Figure()271 fig.add_trace(go.Scatter(x=df['t'], y=df['Substrate'], name='Glucose'))272 fig.add_trace(go.Scatter(x=df['t'], y=df['Biomass'], name='Biomass'))273 fig.add_trace(go.Scatter(x=df['t'], y=df['Acetate'], name='Acetate'))274 fig.update_layout(title=('Simulation of aerobic batch growth of Escherichia coli by acetate cycling'),275 xaxis_title='time (h)',276 yaxis_title='Concentration (g/L)')277 print('print')278 graphJson = json.dumps(fig, cls=plotly.utils.PlotlyJSONEncoder)279 return graphJson280 281 282 283f= Ecoli_Aero()284f.solve()285C= f.solve()[1]286print(C)287 288plt.plot(f.solve()[0], f.solve()[1][:,0])289plt.plot(f.solve()[0], f.solve()[1][:,1])290plt.plot(f.solve()[0], f.solve()[1][:,2])291plt.plot(f.solve()[0], f.solve()[1][:,3])292plt.show()293","Python"
294"Inhibition","simonetacannodelas/BioVL-Library","Rx_Fermentation_Monod-Herbert_Aerobic.py",".py","7444","219","# -*- coding: utf-8 -*-295""""""296Created on Thu Sep 6 13:34:32 2018297 298@author: simoca299""""""300 301 302from scipy.integrate import odeint 303#Package for plotting 304import math 305#Package for the use of vectors and matrix 306import numpy as np307import array as arr308from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas309from matplotlib.figure import Figure310import sys311import os312import matplotlib.pyplot as plt313from matplotlib.ticker import FormatStrFormatter314import glob315from random import sample316import random317import time318 319 320 321class Monod_Herbert:322 def __init__(self, Control=False):323 self.Y_XS = 0.8324 self.Y_OX = 1.05325 self.y_x = 0.5326 self.mu_max = 2.1327 self.Ks = 0.17328 self.kd = 0.21329 self.kla = 1000330 self.O_sat = 0.0755331 332 self.S0 = 18333 self.O0 = 0.0755334 self.X0 = 0.01335 self.V0 = 10336# #parameters for control, default every 1/24 hours:337 self.t_end = 30338 self.t_start = 0339 self.Control = Control340 self.coolingOn = True341 self.steps = (self.t_end - self.t_start)*24342 self.T0 = 30343 self.K_p = 2.31e+01344 self.K_i = 3.03e-01345 self.K_d = -3.58e-03346 self.Tset = 30347 self.u_max = 150348 self.u_min = 0349 def rxn(self, C,t, u):350 #when there is no control, k has no effect351 k=1352 #when cooling is off than u = 0353 if self.coolingOn == False:354 u = 0355 356 if self.Control == True :357 #Cardinal temperature model with inflection: Salvado et al 2011 ""Temperature Adaptation Markedly Determines Evolution within the Genus Saccharomyces""358 #Strain S. cerevisiae PE35 M359 Topt = 30360 Tmax = 45.48361 Tmin = 5.04 362 T = C[5]363 if T < Tmin or T > Tmax:364 k = 0365 else:366 D = (T-Tmax)*(T-Tmin)**2367 E = (Topt-Tmin)*((Topt-Tmin)*(T-Topt)-(Topt-Tmax)*(Topt+Tmin-2*T))368 k = D/E 369 370 #number of components371 n = 3372 m = 3373 #initialize the stoichiometric matrix, s374 s = np.zeros((m,n)) 375 s[0,0] = -1/self.Y_XS376 s[0,1] = -1/self.Y_OX377 s[0,2] = 1378 379 380 s[1,0] = 0381 s[1,1] = 1/self.y_x382 s[1,2] = -1383 384 s[2,0] = 0385 s[2,1] = self.kla386 s[2,2] = 0 387 #initialize the rate vector388 rho = np.zeros((m,1))389 ##initialize the overall conversion vector390 r=np.zeros((n,1))391 rho[0,0] = self.mu_max*(C[0]/(C[0]+self.Ks))*C[2]392 rho[1,0] = self.kd*C[2]393 rho[2,0] = self.kla*(self.O_sat - C[1])394 395 #Developing the matrix, the overall conversion rate is stoichiometric *rates396 r[0,0] = (s[0,0]*rho[0,0])+(s[1,0]*rho[1,0])+(s[2,0]*rho[2,0])397 r[1,0] = (s[0,1]*rho[0,0])+(s[1,1]*rho[1,0])+(s[2,1]*rho[2,0])398 r[2,0] = (s[0,2]*rho[0,0])+(s[1,2]*rho[1,0])+(s[2,2]*rho[2,0])399 400 401 #Solving the mass balances402 dSdt = r[0,0]403 dOdt = r[1,0]404 dXdt = r[2,0]405 dVdt = 0406 if self.Control == True :407 '''408 dHrxn heat produced by cells estimated by yeast heat combustion coeficcient dhc0 = -21.2 kJ/g409 dHrxn = dGdt*V*dhc0(G)-dEdt*V*dhc0(E)-dXdt*V*dhc0(X)410 (when cooling is working) Q = - dHrxn -W ,411 dT = V[L] * 1000 g/L / 4.1868 [J/gK]*dE [kJ]*1000 J/KJ412 dhc0(EtOH) = -1366.8 kJ/gmol/46 g/gmol [KJ/g]413 dhc0(Glc) = -2805 kJ/gmol/180g/gmol [KJ/g]414 415 ''' 416 #Metabolic heat: [W]=[J/s], dhc0 from book ""Bioprocess Engineering Principles"" (Pauline M. Doran) : Appendix Table C.8 417 dHrxndt = dXdt*C[4]*(-21200) #[J/s] + dGdt*C[4]*(15580)- dEdt*C[4]*(29710) 418 #Shaft work 1 W/L1419 W = -1*C[4] #[J/S] negative because exothermic 420 #Cooling just an initial value (constant cooling to see what happens)421 #dQdt = -0.03*C[4]*(-21200) #[J/S] 422 #velocity of cooling water: u [m3/h] -->controlled by PID 423 424 #Mass flow cooling water425 M=u/3600*1000 #[kg/s]426 #Define Tin = 5 C, Tout=TReactor427 #heat capacity water = 4190 J/kgK428 Tin = 5429 #Estimate water at outlet same as Temp in reactor430 Tout = C[5]431 cpc = 4190432 #Calculate Q from Eq 9.47433 Q=-M*cpc*(Tout-Tin) # J/s 434 #Calculate Temperature change435 dTdt = -1*(dHrxndt - Q + W)/(C[4]*1000*4.1868) #[K/s]436 else: 437 dTdt = 0438 return [dSdt, dOdt, dXdt, dVdt, dTdt]439 440 def solve(self):441 #solve normal:442 t = np.linspace(self.t_start, self.t_end, self.steps)443 if self.Control == False :444 u = 0445# fc = 1446 C0 = [self.S0, self.O0, self.X0,self.V0, self.T0]447 C = odeint(self.rxn, C0, t, rtol = 1e-7, mxstep= 500000, args=(u,))448 449 #solve for Control450 else:451 """"""452 PID Temperature Control:453 """"""454 # storage for recording values455 C = np.ones([len(t), 6]) 456 C0 = [self.S0, self.O0, self.X0,self.V0,self.T0]457 C[0] = C0458 self.ctrl_output = np.zeros(len(t)) # controller output459 e = np.zeros(len(t)) # error460 ie = np.zeros(len(t)) # integral of the error461 dpv = np.zeros(len(t)) # derivative of the pv462 P = np.zeros(len(t)) # proportional463 I = np.zeros(len(t)) # integral464 D = np.zeros(len(t)) # derivative465 466 for i in range(len(t)-1):467 #print(t[i])468 #PID control of cooling water469 dt = t[i+1]-t[i]470 #Error471 e[i] = C[i,5] - self.Tset 472 #print(e[i])473 if i >= 1:474 dpv[i] = (C[i,5]-C[i-1,5])/dt475 ie[i] = ie[i-1] + e[i]*dt476 P[i]=self.K_p*e[i]477 I[i]=self.K_i*ie[i]478 D[i]=self.K_d*dpv[i]479 480 self.ctrl_output[i]=P[i]+I[i]+D[i]481 u=self.ctrl_output[i]482 if u>self.u_max:483 u=self.u_max484 ie[i] = ie[i] - e[i]*dt # anti-reset windup485 if u < self.u_min:486 u =self.u_min487 ie[i] = ie[i] - e[i]*dt # anti-reset windup488 #time for solving ODE 489 ts = [t[i],t[i+1]]490 #disturbance491 #if self.t[i] > 5 and self.t[i] < 10:492 # u = 0 493 #solve ODE from last timepoint to new timepoint with old values 494 495 y = odeint(self.rxn, C0, ts, rtol = 1e-7, mxstep= 500000, args=(u,))496 #update C0497 C0 = y[-1]498 #merge y to C499 C[i+1]=y[-1]500 501 return t, C502 503 504 505#f= Monod_Herbert() 506#f.solve() 507#plt.plot(f.solve()[0], f.solve()[1])508#plt.show()509 510 511 512","Python"
513"Inhibition","simonetacannodelas/BioVL-Library","__init__.py",".py","93","8","# -*- coding: utf-8 -*-514""""""515Created on Wed Sep 5 09:51:37 2018516 517@author: bjogut518""""""519 520","Python"
521"Inhibition","simonetacannodelas/BioVL-Library","Sinusoidal oscillating perturbation.py",".py","952","30","#This is the function to add perturbations inside a mechanistic model with concentrations522 523#import the model524 525def disturbances(modelInUse):526 modelInUse.solve()527 t = modelInUse.solve()[0]528 C= modelInUse.solve()[1]529 530 #Initialization of the vectors for the superposition of the perturbations531 PV =[]532 PV2 = []533 PV3 = []534 535 #Creation of the vector to collect the data with the noise/perturbation536 C_noise = np.zeros((C.shape))537 import numpy as np538 539 #sinusoidal oscillation of each time - Dependance of time540 for i in range(len(t)):541 PV.append((((math.sin(t[i]*83))/47)+(math.cos(t[i]*11))/23))542 PV2.append(0.6*(((math.sin(t[i] * 61)/31) + (math.cos(t[i] * 43)) / 61)))543 PV3.append(0.33*(((math.sin(t[i] * 41))/19) + (math.cos(t[i] * 61)) / 89))544 545 #stablishing the perturbation inside the concentration546 for i in range(len(C[0])):547 C_noise [:, i] = C[:, i] + C[:, i]* PV + C[:, i]*PV2 + C[:, i]*PV3548 549 return C_noise550","Python"
551"Inhibition","simonetacannodelas/BioVL-Library","Rx_Fermentation_Scerevisiae_Glucose_Aerobic_Batch.py",".py","8529","234","# -*- coding: utf-8 -*-552""""""553Created on Thu Sep 6 13:34:32 2018554 555@author: bjogut and simoca556""""""557 558 559from scipy.integrate import odeint 560#Package for plotting 561import math 562#Package for the use of vectors and matrix 563import numpy as np564import array as arr565from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas566from matplotlib.figure import Figure567import sys568import os569import matplotlib.pyplot as plt570from matplotlib.ticker import FormatStrFormatter571import glob572from random import sample573import random574import time575 576 577class SCerevisiae_Aero:578 def __init__(self, Control = False):579 self.Yox_XG = 0.8580 self.Yred_XG = 0.05581 self.Yox_XE = 0.72582 self.Y_OG = 1.067583 self.Y_EG = 0.5584 self.Y_OE = 1.5585 self.q_g = 3.5586 self.q_o = 0.37587 self.q_e = 0.32588 self.t_lag = 4.66589 self.Kg = 0.17590 self.Ke = 0.56591 self.Ko = 0.0001592 self.Ki = 0.31593 self.O_sat = 0.00755594 self.kla = 1004595 596 self.G0 = 18597 self.E0 = 0.0598 self.O0 = 0.00755599 self.X0 = 0.1600 601 self.t_end = 30602 self.t_start = 0603 self.V0 = 2604 605 #parameters for control, default every 1/24 hours:606 self.Control = Control607 self.coolingOn = True608 self.Contamination=False609 self.steps = (self.t_end - self.t_start)*24610 self.T0 = 30611 self.K_p = 2.31e+01612 self.K_i = 3.03e-01613 self.K_d = -3.58e-03614 self.Tset = 30615 self.u_max = 150616 self.u_min = 0617 618 def rxn(self, C,t , u, fc):619 #when there is no control, k has no effect620 k=1621 #when cooling is off than u = 0622 if self.coolingOn == False:623 u = 0624 if self.Contamination == True:625 fc=np.random.randint(0,10)626 fc=fc/17627 628 if self.Control == True :629 #Cardinal temperature model with inflection: Salvado et al 2011 ""Temperature Adaptation Markedly Determines Evolution within the Genus Saccharomyces""630 #Strain S. cerevisiae PE35 M631 Topt = 30632 Tmax = 45.48633 Tmin = 5.04 634 T = C[5]635 if T < Tmin or T > Tmax:636 k = 0637 else:638 D = (T-Tmax)*(T-Tmin)**2639 E = (Topt-Tmin)*((Topt-Tmin)*(T-Topt)-(Topt-Tmax)*(Topt+Tmin-2*T))640 k = D/E 641 642 #number of components643 n = 4644 m = 4645 #initialize the stoichiometric matrix, s646 s = np.zeros((m,n)) 647 s[0,0] = -1648 s[0,1] = 0649 s[0,2] = -self.Y_OG650 s[0,3] = self.Yox_XG651 652 s[1,0] = -1653 s[1,1] = self.Y_EG654 s[1,2] = 0655 s[1,3] = self.Yred_XG656 657 s[2,0] = 0658 s[2,1] = -1659 s[2,2] = -self.Y_OE660 s[2,3] = self.Yox_XE661 662 s[3,0] = 0663 s[3,1] = 0664 s[3,2] = 1665 s[3,3] = 0666 #initialize the rate vector667 rho = np.zeros((4,1))668 ##initialize the overall conversion vector669 r=np.zeros((4,1))670 rho[0,0] = k*((1/self.Y_OG)*min(self.q_o*(C[2]/(C[2]+self.Ko)),self.Y_OG*(self.q_g*(C[0]/(C[0]+self.Kg)))))*C[3]671 rho[1,0] = k*((1-math.exp(-t/self.t_lag))*((self.q_g*(C[0]/(C[0]+self.Kg)))-(1/self.Y_OG)*min(self.q_o*(C[2]/(C[2]+self.Ko)),self.Y_OG*(self.q_g*(C[0]/(C[0]+self.Kg))))))*C[3]672 rho[2,0] = k*((1/self.Y_OE)*min(self.q_o*(C[2]/(C[2]+self.Ko))-(1/self.Y_OG)*min(self.q_o*(C[2]/(C[2]+self.Ko)),self.Y_OG*(self.q_g*(C[0]/(C[0]+self.Kg)))),self.Y_OE*(self.q_e*(C[1]/(C[1]+self.Ke))*(self.Ki/(C[0]+self.Ki)))))*C[3]673 rho[3,0] = self.kla*(self.O_sat - C[2])674 675 #Developing the matrix, the overall conversion rate is stoichiometric *rates676 r[0,0] = (s[0,0]*rho[0,0])+(s[1,0]*rho[1,0])+(s[2,0]*rho[2,0])+(s[3,0]*rho[3,0])677 r[1,0] = (s[0,1]*rho[0,0])+(s[1,1]*rho[1,0])+(s[2,1]*rho[2,0])+(s[3,1]*rho[3,0])678 r[2,0] = (s[0,2]*rho[0,0])+(s[1,2]*rho[1,0])+(s[2,2]*rho[2,0])+(s[3,2]*rho[3,0])679 r[3,0] = (s[0,3]*rho[0,0])+(s[1,3]*rho[1,0])+(s[2,3]*rho[2,0])+(s[3,3]*rho[3,0])680 681 #Solving the mass balances682 dGdt = r[0,0]683 dEdt = r[1,0]*fc684 dOdt = r[2,0]685 dXdt = r[3,0]686 dVdt = 0687 if self.Control == True :688 '''689 dHrxn heat produced by cells estimated by yeast heat combustion coeficcient dhc0 = -21.2 kJ/g690 dHrxn = dGdt*V*dhc0(G)-dEdt*V*dhc0(E)-dXdt*V*dhc0(X)691 (when cooling is working) Q = - dHrxn -W ,692 dT = V[L] * 1000 g/L / 4.1868 [J/gK]*dE [kJ]*1000 J/KJ693 dhc0(EtOH) = -1366.8 kJ/gmol/46 g/gmol [KJ/g]694 dhc0(Glc) = -2805 kJ/gmol/180g/gmol [KJ/g]695 696 ''' 697 #Metabolic heat: [W]=[J/s], dhc0 from book ""Bioprocess Engineering Principles"" (Pauline M. Doran) : Appendix Table C.8 698 dHrxndt = dXdt*C[4]*(-21200) #[J/s] + dGdt*C[4]*(15580)- dEdt*C[4]*(29710) 699 #Shaft work 1 W/L1700 W = 1*C[4] #[J/S] negative because exothermic 701 #Cooling just an initial value (constant cooling to see what happens)702 #dQdt = -0.03*C[4]*(-21200) #[J/S] 703 #velocity of cooling water: u [m3/h] -->controlled by PID 704 705 #Mass flow cooling water706 M=u/3600*1000 #[kg/s]707 #Define Tin = 5 C, Tout=TReactor708 #heat capacity water = 4190 J/kgK709 Tin = 5710 #Estimate water at outlet same as Temp in reactor711 Tout = C[5]712 cpc = 4190713 #Calculate Q from Eq 9.47714 Q=-M*cpc*(Tout-Tin) # J/s 715 #Calculate Temperature change716 dTdt = -1*(dHrxndt - Q + W)/(C[4]*1000*4.1868) #[K/s]717 else: 718 dTdt = 0719 return [dGdt,dEdt,dOdt,dXdt,dVdt, dTdt]720 721 def solve(self):722 #solve normal:723 t = np.linspace(self.t_start, self.t_end, self.steps)724 if self.Control == False :725 u = 0726 fc= 1727 C0 = [self.G0, self.E0, self.O0, self.X0,self.V0,self.T0]728 C = odeint(self.rxn, C0, t, rtol = 1e-7, mxstep= 500000, args=(u,fc,))729 730 #solve for Control731 else:732 fc=0733 """"""734 PID Temperature Control:735 """"""736 # storage for recording values737 C = np.ones([len(t), 6]) 738 C0 = [self.G0, self.E0, self.O0, self.X0,self.V0,self.T0]739 C[0] = C0740 self.ctrl_output = np.zeros(len(t)) # controller output741 e = np.zeros(len(t)) # error742 ie = np.zeros(len(t)) # integral of the error743 dpv = np.zeros(len(t)) # derivative of the pv744 P = np.zeros(len(t)) # proportional745 I = np.zeros(len(t)) # integral746 D = np.zeros(len(t)) # derivative747 748 for i in range(len(t)-1):749 #print(t[i])750 #PID control of cooling water751 dt = t[i+1]-t[i]752 #Error753 e[i] = C[i,5] - self.Tset 754 #print(e[i])755 if i >= 1:756 dpv[i] = (C[i,5]-C[i-1,5])/dt757 ie[i] = ie[i-1] + e[i]*dt758 P[i]=self.K_p*e[i]759 I[i]=self.K_i*ie[i]760 D[i]=self.K_d*dpv[i]761 762 self.ctrl_output[i]=P[i]+I[i]+D[i]763 u=self.ctrl_output[i]764 if u>self.u_max:765 u=self.u_max766 ie[i] = ie[i] - e[i]*dt # anti-reset windup767 if u < self.u_min:768 u =self.u_min769 ie[i] = ie[i] - e[i]*dt # anti-reset windup770 #time for solving ODE 771 ts = [t[i],t[i+1]]772 #disturbance773 #if self.t[i] > 5 and self.t[i] < 10:774 # u = 0 775 #solve ODE from last timepoint to new timepoint with old values 776 777 y = odeint(self.rxn, C0, ts, rtol = 1e-7, mxstep= 500000, args=(u,fc,))778 #update C0779 C0 = y[-1]780 #merge y to C781 C[i+1]=y[-1]782 return t, C783 784","Python"
785"Inhibition","simonetacannodelas/BioVL-Library","Rx_Fermentation_Monod-Herbert_Anaerobic.py",".py","8180","233","# -*- coding: utf-8 -*-786""""""787Created on Thu Sep 6 13:34:32 2018788 789@author: simoca790""""""791 792 793from scipy.integrate import odeint 794#Package for plotting 795import math 796#Package for the use of vectors and matrix 797import numpy as np798import array as arr799from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas800from matplotlib.figure import Figure801import sys802import os803import matplotlib.pyplot as plt804from matplotlib.ticker import FormatStrFormatter805import glob806from random import sample807import random808import time809 810 811class Monod_Herbert_Anaero:812 def __init__(self, Control=False):813 self.Y_XS = 0.8814 self.Y_OX = 1.05815 self.Y_PX = 2816 self.y_x = 0.5817 self.mu_max = 2.1818 self.Ks = 0.17819 self.kd = 0.21820 self.kla = 1000821 self.O_sat = 0.0755822 823 self.S0 = 18824 self.X0 = 0.1825 self.V0 = 10826 self.P0=0827# #parameters for control, default every 1/24 hours:828 self.t_end = 30829 self.t_start = 0830 self.Control = Control831 self.coolingOn = True832 self.steps = (self.t_end - self.t_start)*24833 self.T0 = 30834 self.K_p = 2.31e+01835 self.K_i = 3.03e-01836 self.K_d = -3.58e-03837 self.Tset = 30838 self.u_max = 150839 self.u_min = 0840# def compounds(self):841# self.compound = ['Substrate','Product','Biomass']842# return self.compound843 def rxn(self, C,t, u):844 #when there is no control, k has no effect845 k=1846 #when cooling is off than u = 0847 if self.coolingOn == False:848 u = 0849 850 if self.Control == True :851 #Cardinal temperature model with inflection: Salvado et al 2011 ""Temperature Adaptation Markedly Determines Evolution within the Genus Saccharomyces""852 #Strain S. cerevisiae PE35 M853 Topt = 30854 Tmax = 45.48855 Tmin = 5.04 856 T = C[5]857 if T < Tmin or T > Tmax:858 k = 0859 else:860 D = (T-Tmax)*(T-Tmin)**2861 E = (Topt-Tmin)*((Topt-Tmin)*(T-Topt)-(Topt-Tmax)*(Topt+Tmin-2*T))862 k = D/E 863 864 #number of components865 self.s = np.zeros((2,3))866 self.rho=np.zeros((2,1)) 867 self.s[0,2]=1868 self.s[0,0]=(-1/self.Y_XS)869 self.s[0,1] = (1/self.Y_PX)870 self.s[1,2]=-1871 872 self.rho[0,0]=((self.mu_max*C[0])/(C[0]+self.Ks))*C[2]873 self.rho[1,0]=self.kd*C[2]874 875# print(self.rho)876# print(self.s)877 self.r= np.zeros((3,1))878 self.r[0,0]= self.s[0,0]*self.rho[0,0]+self.s[1,0]*self.rho[1,0]879 self.r[1,0]= self.s[0,1]*self.rho[0,0]+self.s[1,1]*self.rho[1,0]880 self.r[2,0]= self.s[0,2]*self.rho[0,0]+self.s[1,2]*self.rho[1,0]881 dSdt = self.r[0,0]882 dPdt = self.r[1,0]883 dXdt = self.r[2,0]884 dVdt = 0885# 886# n = 3887# m = 3888# #initialize the stoichiometric matrix, s889# s = np.zeros((m,n)) 890# s[0,0] = -1/self.Y_XS891# s[0,1] = -1/self.Y_OX892# s[0,2] = 1893# 894# 895# s[1,0] = 0896# s[1,1] = 1/self.y_x897# s[1,2] = -1898# 899# s[2,0] = 0900# s[2,1] = self.kla901# s[2,2] = 0 902# #initialize the rate vector903# rho = np.zeros((m,1))904# ##initialize the overall conversion vector905# r=np.zeros((n,1))906# rho[0,0] = self.mu_max*(C[0]/(C[0]+self.Ks))*C[2]907# rho[1,0] = self.kd*C[2]908# rho[2,0] = self.kla*(self.O_sat - C[1])909# 910# #Developing the matrix, the overall conversion rate is stoichiometric *rates911# r[0,0] = (s[0,0]*rho[0,0])+(s[1,0]*rho[1,0])+(s[2,0]*rho[2,0])912# r[1,0] = (s[0,1]*rho[0,0])+(s[1,1]*rho[1,0])+(s[2,1]*rho[2,0])913# r[2,0] = (s[0,2]*rho[0,0])+(s[1,2]*rho[1,0])+(s[2,2]*rho[2,0])914# 915#916# #Solving the mass balances917# dSdt = r[0,0]918# dOdt = r[1,0]919# dXdt = r[2,0]920# dVdt = 0921 if self.Control == True :922 '''923 dHrxn heat produced by cells estimated by yeast heat combustion coeficcient dhc0 = -21.2 kJ/g924 dHrxn = dGdt*V*dhc0(G)-dEdt*V*dhc0(E)-dXdt*V*dhc0(X)925 (when cooling is working) Q = - dHrxn -W ,926 dT = V[L] * 1000 g/L / 4.1868 [J/gK]*dE [kJ]*1000 J/KJ927 dhc0(EtOH) = -1366.8 kJ/gmol/46 g/gmol [KJ/g]928 dhc0(Glc) = -2805 kJ/gmol/180g/gmol [KJ/g]929 930 ''' 931 #Metabolic heat: [W]=[J/s], dhc0 from book ""Bioprocess Engineering Principles"" (Pauline M. Doran) : Appendix Table C.8 932 dHrxndt = dXdt*C[4]*(-21200) #[J/s] + dGdt*C[4]*(15580)- dEdt*C[4]*(29710) 933 #Shaft work 1 W/L1934 W = -1*C[4] #[J/S] negative because exothermic 935 #Cooling just an initial value (constant cooling to see what happens)936 #dQdt = -0.03*C[4]*(-21200) #[J/S] 937 #velocity of cooling water: u [m3/h] -->controlled by PID 938 939 #Mass flow cooling water940 M=u/3600*1000 #[kg/s]941 #Define Tin = 5 C, Tout=TReactor942 #heat capacity water = 4190 J/kgK943 Tin = 5944 #Estimate water at outlet same as Temp in reactor945 Tout = C[5]946 cpc = 4190947 #Calculate Q from Eq 9.47948 Q=-M*cpc*(Tout-Tin) # J/s 949 #Calculate Temperature change950 dTdt = -1*(dHrxndt - Q + W)/(C[4]*1000*4.1868) #[K/s]951 else: 952 dTdt = 0953 return [dSdt, dPdt, dXdt, dVdt, dTdt] 954 955 def solve(self):956 #solve normal:957 t = np.linspace(self.t_start, self.t_end, self.steps)958 if self.Control == False :959 u = 0960 C0 = [self.S0, self.P0, self.X0,self.V0, self.T0]961 C = odeint(self.rxn, C0, t, rtol = 1e-7, mxstep= 500000, args=(u,))962 963 #solve for Control964 else:965 """"""966 PID Temperature Control:967 """"""968 # storage for recording values969 C = np.ones([len(t), 6]) 970 C0 = [self.S0, self.P0, self.X0,self.V0,self.T0]971 self.ctrl_output = np.zeros(len(t)) # controller output972 e = np.zeros(len(t)) # error973 ie = np.zeros(len(t)) # integral of the error974 dpv = np.zeros(len(t)) # derivative of the pv975 P = np.zeros(len(t)) # proportional976 I = np.zeros(len(t)) # integral977 D = np.zeros(len(t)) # derivative978 979 for i in range(len(t)-1):980 #print(t[i])981 #PID control of cooling water982 dt = t[i+1]-t[i]983 #Error984 e[i] = C[i,5] - self.Tset 985 #print(e[i])986 if i >= 1:987 dpv[i] = (C[i,5]-C[i-1,5])/dt988 ie[i] = ie[i-1] + e[i]*dt989 P[i]=self.K_p*e[i]990 I[i]=self.K_i*ie[i]991 D[i]=self.K_d*dpv[i]992 993 self.ctrl_output[i]=P[i]+I[i]+D[i]994 u=self.ctrl_output[i]995 if u>self.u_max:996 u=self.u_max997 ie[i] = ie[i] - e[i]*dt # anti-reset windup998 if u < self.u_min:999 u =self.u_min1000 ie[i] = ie[i] - e[i]*dt # anti-reset windup1001 #time for solving ODE 1002 ts = [t[i],t[i+1]]1003 #disturbance1004 #if self.t[i] > 5 and self.t[i] < 10:1005 # u = 0 1006 #solve ODE from last timepoint to new timepoint with old values 1007 1008 y = odeint(self.rxn, C0, ts, rtol = 1e-7, mxstep= 500000, args=(u,))1009 #update C01010 C0 = y[-1]1011 #merge y to C1012 C[i+1]=y[-1]1013 1014 return t, C1015 1016 1017","Python"
1018"Inhibition","simonetacannodelas/BioVL-Library","Rx_Fermentation_Ecoli_Glucose_Aerobic_Batch.py",".py","10981","296","# -*- coding: utf-8 -*-1019""""""1020Created on Friday Jan 24 13:34:32 20201021 1022@author: simoca1023""""""1024from scipy.integrate import odeint 1025#Package for plotting 1026import math 1027#Package for the use of vectors and matrix 1028import numpy as np1029import pandas as pd1030import array as arr1031from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas1032from matplotlib.figure import Figure1033import sys1034import os1035import matplotlib.pyplot as plt1036from matplotlib.ticker import FormatStrFormatter1037import glob1038from random import sample1039import random1040import time1041import plotly1042import plotly.graph_objs as go1043import json1044from plotly.subplots import make_subplots1045 1046 1047class Ecoli_Aero:1048 def __init__(self, Control = False):1049 self.Kap=0.5088 #g/L1050 self.Ksa=0.01281051 self.Kia=1.2602 #g/L1052 self.Ks= 0.0381 #g/L1053 self.Kis = 1.8383 #g/L affinity constant1054 self.Ko = 0.00011055 self.qAcmax= 0.1148 #g/gh1056 self.qm=0.0133 #g/gh1057 self.qOmax= 13.4*31.9988/1000 #g/gh1058 self.qSmax= 0.635 #g/gh1059 self.Yas= 0.8938 #g/g1060 self.Yoa= 0.5221 #g/g1061 self.Yos= 1.5722 #g/g1062 self.Yxa= 0.5794 #g/g1063 self.Yem = 0.5321 #g/g1064 self.Yxsof = 0.229 # g/g1065 self.pAmax= 0.2286 #gA/gXh1066 self.kla= 2201067 self.H= 1400 #Henry constant1068 1069 self.G0 = 4.941070 self.A0 = 0.01291071 self.O0 = 0.09771072 self.X0 = 0.171073 self.tau= 35 #response time1074 self.t_start = 01075 self.V0 = 21076 self.F0 = 01077 self.SFR = 1.51078 self.t_expfb_start = 01079 self.t_constfb_start = 1.51080 self.t_end = 301081 1082 #parameters for control, default every 1/24 hours:1083 self.Control = Control1084 self.coolingOn = True1085 self.Contamination=False1086 self.steps = (self.t_end - self.t_start)*241087 self.T0 = 351088 self.K_p = 2.31e+011089 self.K_i = 11090 self.K_d = 01091 self.Tset = 301092 self.u_max = 1501093 self.u_min = 01094 1095 def update_param_value(self, param, new_value):1096 1097 class_name = str(self.__class__).split('.')[1].split(""'"")[0]1098 param_current_value = self.__dict__.get(param, None)1099 1100 if param_current_value is None:1101 print(""Class {} does not contain param with name {}"".format(class_name, param))1102 return1103 1104 self.__dict__[param] = new_value1105 print(""Value of param {} in class {} updated to {}"".format(param, class_name, new_value))1106 1107 1108 def rxn(self, C,t , u, fc):1109 #when there is no control, k has no effect1110 k=11111 #when cooling is off than u = 01112 if self.coolingOn == False:1113 u = 01114 if self.Contamination == True:1115 fc=np.random.randint(0,10)1116 fc=fc/171117 1118 if self.Control == True :1119 #Cardinal temperature model with inflection: Salvado et al 2011 ""Temperature Adaptation Markedly Determines Evolution within the Genus Saccharomyces""1120 #Strain E.coli W3101121 Topt = 351122 Tmax = 45.481123 Tmin = 101124 T = C[5]1125 if T < Tmin or T > Tmax:1126 k = 01127 else:1128 D = (T-Tmax)*(T-Tmin)**21129 E = (Topt-Tmin)*((Topt-Tmin)*(T-Topt)-(Topt-Tmax)*(Topt+Tmin-2*T))1130 k = D/E1131 # Volume balance1132 if (t >= self.t_expfb_start):1133 if (t < self.t_constfb_start):1134 Fin = self.F0 * math.exp(self.SFR * (t - self.t_expfb_start))1135 Fout = 01136 else:1137 Fin = self.F0 * math.exp(self.SFR * (self.t_constfb_start - self.t_expfb_start))1138 Fout = 01139 else:1140 Fin = 01141 Fout = 01142 1143 F = Fin - Fout1144 1145 qS = (self.qSmax/(1+C[1]/self.Kia))*(C[0]/(C[0]+self.Ks))1146 qSof = self.pAmax*(qS/(qS+self.Kap))1147 pA = qSof*self.Yas1148 qSox = (qS-qSof)*(C[2]/(C[2]+self.Ko))1149 qSan = (qSox-self.qm)*self.Yem*(0.488/0.391)1150 qsA = (self.qAcmax/(1+(qS/self.Kis)))*(C[1]/(C[1]+self.Ksa))1151 qA = pA - qsA1152 mu = (qSox - self.qm)*self.Yem+ qsA*self.Yxa + qSof*self.Yxsof1153 qO = self.Yos*(qSox-qSan)+qsA*self.Yoa1154 1155 #Solving the mass balances1156 dGdt = -(qS*C[3])1157 dAdt = qA*C[3]1158 dOdt = self.kla*(self.O0-C[2])1159 dXdt = (mu)*C[3]1160 dVdt = 01161 if self.Control == True :1162 '''1163 dHrxn heat produced by cells estimated by yeast heat combustion coeficcient dhc0 = -21.2 kJ/g1164 dHrxn = dGdt*V*dhc0(G)-dEdt*V*dhc0(E)-dXdt*V*dhc0(X)1165 (when cooling is working) Q = - dHrxn -W ,1166 dT = V[L] * 1000 g/L / 4.1868 [J/gK]*dE [kJ]*1000 J/KJ1167 dhc0(EtOH) = -1366.8 kJ/gmol/46 g/gmol [KJ/g]1168 dhc0(Glc) = -2805 kJ/gmol/180g/gmol [KJ/g]1169 1170 ''' 1171 #Metabolic heat: [W]=[J/s], dhc0 from book ""Bioprocess Engineering Principles"" (Pauline M. Doran) : Appendix Table C.8 1172 dHrxndt = dXdt*C[4]*(-21200) #[J/s] + dGdt*C[4]*(15580)- dEdt*C[4]*(29710) 1173 #Shaft work 1 W/L11174 W = 1*C[4] #[J/S] negative because exothermic 1175 #Cooling just an initial value (constant cooling to see what happens)1176 #dQdt = -0.03*C[4]*(-21200) #[J/S] 1177 #velocity of cooling water: u [m3/h] -->controlled by PID 1178 1179 #Mass flow cooling water1180 M=u/3600*1000 #[kg/s]1181 #Define Tin = 5 C, Tout=TReactor1182 #heat capacity water = 4190 J/kgK1183 Tin = 51184 #Estimate water at outlet same as Temp in reactor1185 Tout = C[5]1186 cpc = 41901187 #Calculate Q from Eq 9.471188 Q=-M*cpc*(Tout-Tin) # J/s 1189 #Calculate Temperature change1190 dTdt = -1*(dHrxndt - Q + W)/(C[4]*1000*4.1868) #[K/s]1191 else: 1192 dTdt = 01193 return [dGdt,dAdt,dOdt,dXdt,dVdt, dTdt]1194 1195 def solve(self):1196 #solve normal:1197 t = np.linspace(self.t_start, self.t_end, self.steps)1198 if self.Control == False :1199 u = 01200 fc= 1