Metabolic Engineering
Prof. Amit Ghosh
School of Energy Science and Engineering
Indian Institute of Technology, Kharagpur
Lecture - 20
Optimization in MATLAB
Welcome to metabolic engineering course. Today we will learn about the optimization technique in MATLAB. Some basics in MATLAB are required we can handle matrix.
(Refer Slide Time: 00:41)
Today, we will cover these topics the matrix operation, especially the matrix operation where you can do a lot of operation in matrix that will go into that we are going to learn and then followed by network optimization and more practical example while we will consider the E. coli core metabolic network optimization. How you can optimize the entire network while you have several more reaction.
(Refer Slide Time: 01:15)
So, we will start with the matrix multiplication like how we can write a matrix in a MATLAB matrix A is given over here and then matrix B how we can multiply matrix A and matrix B, that we are going to align and also B multiplied by A and then how you can calculate the null space of A. Let us start with let us open the MATLAB in your computer and we will start doing some operation here.
(Refer Slide Time: 01:53)
So, for example, let us define a matrix how you can define a matrix. Suppose I define the matrix A which I have already told that matrix A is equal to third bracket and then 1 2 1 2 then semicolon and 3 4 -3 4 1 3 1 3 then again semicolon. So, every row you are putting at the end you are putting a semicolon and then the last row 2 4 2 4 and then semicolon enter. Now, you can see the matrix we have entered in the MATLAB. So, exactly where the matrix A which I have defined that is 1 2 1 2 3 4 -3 4 and then 1 3 1 3 2 4 2 4.
So, this matrix you can print the row. Suppose I want to print the size of the matrix, what is the size of the matrix the size of the matrix you can just put size in the bracket you put A matrix A. So, your size is a 4 by 4 matrix, the dimension is M = 4 and N = 4, that is 4 rows and 4 columns and then I want to print the contents of the; suppose I want to print this second row.
So, second row if I want to print then it is just A then for then since I want to print this second row 2, then the colon and then it will print the sub array that is the second row 3 4 -3 4. Then if I want to print the submatrix here like if I were if I just want to print this one this sub array this starting from 4, -3, 4, 3, 1, 3, 4, 2, 4 this sub matrix we start from 4 and then and the last element is 4. So, how you can actually print the subarray that also you can actually print, or you can store that is A bracket 1, 2, 2 and then 2 2 4.
So, it is actually printing since I consider first row second element 1 1 2 1 2 2 that is row 1 element 2 that then that is 2 1 2 and then we have 2 2 4 then it is from 4 2 -4 -3 4. So, this will be this way you can print different sub array from the matrix and if you enter the second matrix which I want to multiply the second matrix is that B to the if I enter the matrix B and then B will be 1, 2, 0, 0 semicolon 4, 0, 1, 3 and then 3, 0, 0, 3; 1, 2, 1, 2.
Then you are closing the matrix, so, this is your second B matrix. So, how do you multiply matrix A and B it is just A then we have the star multiply sign from your keyboard and then B we just multiply A and B and your new matrix will be just A multiplied by B. So, your new matrix is 14, 6, 4, 13, 14, 14, 8, 11, 19, 8, 6, 18, 28, 12, 8, 20. So, you can check in your MATLAB you can multiply A and B very easily suppose, I multiply B and B with A, B multiplied A then enter now, you can see that your matrix is different 7, 7, 10, -5, 10.
So, A multiplied B is not equal to B multiplied A, A and B are not same. So, that is why they are different and now, I want to calculate the find the basis of the null space. So, how we can calculate the null space of a matrix is basically null A. So, now we can see the null space vector is basically 0, 0.7071, 0.000, – 0.7071. So, this is your null space basis vector. So, using that basis vector you can draw your solution space lies in under this basis vector.
Similarly, if I if you type null B then it will give you the null space basis vector for this matrix B is 0 this 4 by 0 is a basically 0 matrix all the elements are 0. This way you can
calculate the basis vector for different matrix just by typing null, within bracket you can specify the matrix. So, now, today we learn about how you can actually enter a matrix A and also matrix B and how to multiply A and B and also B and A and we have seen that AB is not equal to BA. The multiplication is not commutative and then we also calculate the basis for null space A.
(Refer Slide Time: 09:46)
Now, we want to say solve another problem that is a linear equation how we can solve this linear equation in MATLAB. So, this linear equation is very similar to fluxes that you have dx / dt equal to those linear equation can be represented by a matrix form. So, all the linear equation you get from s dot v = 0 can be used to form a can I can do matrix algebra using those equations.
So, this is one example where you see the matrix2 0 3 these are these tachymetric coefficients and the fluxes are maybe x1, x2, x3 and generally for the steady state we are assumed that the on the right-hand side is 6, -1, -16 when it is in a steady state it is 0. So, now, how we can solve this equation how we can get the value of x1, x2, x3 so, you can do it by hand or using MATLAB there is a trick you can calculate in one very easily however, you can calculate the value of x1, x2, x3 very easily how to do that.
(Refer Slide Time: 11:05)
So, this is A this linear set of equations and then I define the left-hand matrix as A I can be defined as 2 0 3 and then 1 4 -1 and -5, 1, 1 this is matrix on the left-hand side 20314 -1 -511. Now, we want to define another matrix B, B = 6 -1 -16. So, this is your next matrix is the column 1 column matrix it is 3 by 1 now, we want to define another vector that is X. X = x1 semicolon x2, x3 this is your variable.
So, since it is a variable in a MATLAB we have to actually put under quotation, so, these x1, x2, x3 this is another added variable x1, x2, x3 now you want to say solve this linear equation, so, you can calculate X as X = B B/A, these are the value for x1 x2 x3, your X will be X is a column matrix. So, to solve X you have to do A/B. So, your AX x3 is 0, x2 is -1, x1 is 3. So, let us check. So, your first equation is 2x1 + 3x3= 6.
This is your first equation and then you can solve it by hand also you can solve it by hand and so, 2x1 + 3x3 = 6, so, just substitute that value, so, x1 is 3 the 3 into 6 and then 3 into 0. So, it is equal to 0. Correct. So, 6 = 6. So, first equation is correct, this is the second equation is x1 - x1 + 4x2 so, matrix are multiplied in this fashion multiply this and with this, so, it can actually A multiplied into X. So now your next equation will be x1 + 4x2 - x3.
So, let us have a substitute the value your x1 is 3, x2 is -1and x3 is 0. So, that is equal to -1. So, if we substitute the value of x1 x2 x3 which are given here, I show your x1 is 3 x2 is -1. So, 4 into -1 -4x3 is 0. So, it is basically it is say 4 - x1 – 4x2 - x3 that is equal to -1. So, it is matching. So, x1 is 3 -4 into x2 that is 1. So, 1 -4 3 - 4 it is becoming -1. So, it is correct is also matching.
Your first equation is also my this is the first equation you are getting by multiplying this row with this column and then this row with this column is -5x1 + x2 + x3 = -16. So, you are multiplying this row with this column. So, this way, so, you can see that the x1 is 3, so 3 into 5 15 -15 + x2 -1 that is 15-1and x3 is 0. So, it will be -15. So, is the value solutions are actually following this equation.
So, this way, you can actually multiply and the formula is, so if you want to, if this is so matrix operation is that you multiply this row with this column, this row with this column with this row with this column. So, basically, your first equation will become 2x1 + 3x3 = 6 and your second equation is basically x1 + 4x2 -x3 = -1 and then your third equation is basically 5x1 + x2 + x3 = -16.
So, these are the 3 linear equation you are getting from this equation and from the MATLAB what we saw that if you define these as matrix B matrix A and this is X and this is B. Then the formula that you use for calculating X, X is basically A slash B. So, this is the formula that we use in MATLAB. So, this is the matrix we want to calculate. So, they were the MATLAB for formula that we used basically A this matrix the slash, the symbol you should remember this is this slash, which you can find in the keyboard.
And then then B, B is basically the numerical value that is 6 -1 -16. So, this way, you can solve any linear equation in a matrix form. So, this using MATLAB, you can solve any set of linear equation, which is x1 x2 x3 and before and I have already shown that what are the value of x1 and x2.
(Refer Slide Time: 20:41)
So from your calculation, that is A, B, A / B, so these becoming your solution. So x1 and the solution is basically given by 3-1 0. So, your solution for x1 is x1 = 3 x2 = -1 and x3 = 0. So, these are the 3 values you are getting just by using this command. So, if you put this command A / B then you get this value x3 x1 = 3, x2 = -1 and x3 = 0. So, if you substitute this value, the first equation, x1 = 3, 6 and x3 = 0, so 6 = 6.
So, first equation is valid, then you put x1 = 3 x2 = -1, -1 into 4, that is -4, so 3 -4 -1. x3 is anyway 0, all x3 are 0, then the third equation, x1, if you multiply 5 with 3 15 -15, -1, so it becomes -16. So, all 3 equations are actually true, if you have a solution, x1 = 3, x2 = -1 and x3 = 0. So, this way, you can easily solve any linear set of equations that you get and then try to get a solution. Now we will learn how to actually calculate a determinant of a matrix.
So, matrix determinant you can calculate determinant of a matrix, this is a simple matrix1 3 2 3 what is the value here, so it is 1 3 2 3. So, you define a matrix1 3 2 3. So, this matrix determinant you can calculate det within the bracket A. So, the determinant of A is -3, let us calculate by hand with A -3 is correct or wrong.
(Refer Slide Time: 23:11)
So, the determinant of a matrix is basically so 1 multiplied by 3 that is -2 into 3, 3-6, so your determinant is -3. So, if from my hand also, we saw that the determinant is -3 when using MATLAB also we saw, the determinant is -3. So, fine so you can see how to multiply 1 into 3, this 1 and this 1 first, it is coming here and then 2 into 3, so -2 into this way you can calculate and the determinant of matrix.
(Refer Slide Time: 24:06)
Now we will learn about how to calculate the eigenvalue and eigenvector. So eigenvalues of a given matrix, say this one that is G equal to 2 7 7 2 and how we can calculate the eigenvalues and eigenvectors. So, using by hand also you can calculate but in MATLAB it is very easy. Just to enter one command in the command is eigen eig then the matrix you define. So, this is the command you are going to use. So, V, L. So, V is a eigen vector and L is a eigen value, I using this command, you are going to get the first you enter the matrix and then you calculate the vector and the eigen values.
(Refer Slide Time: 25:10)
So, let us do that we will 2 7, G = 2 7 7 2. Now, you just use that command that is V L close the bracket and then it equal to eig within bracket G. So, now you can see the eigenvalues are -5 and 9 eigenvalues are always stored in a diagonal. So, these are the eigenvector, eigenvalue is -5 9 and then corresponding eigen vectors are -1, -1 you can the constant term you can put outside will be a -1 1 1 1.
So, finally, what do you get the eigenvalues. So, your lambda 1 through n, you solve by hand your lambda 1 will be -1 and lambda 2 = 6 -5 and 9 lambda 1 is -5, lambda 2 is 9 and corresponding eigenvector are 1 1 and you just find 7 0 7 you can bring outside so your eigen vectors are x1 and x2 given by 1 -1. So, this way can calculate the eigenvalues and eigenvectors. Now we will learn about the optimization scheme.
(Refer Slide Time: 27:25)
Now we will start with the network metabolic network optimization. So, we learned about metabolic network today we are going to do this calculation how to optimize the network. So, earlier we saw that how we can construct the network. This is the network while we have the metabolic which is coming in and then there is a reaction v 1 v 2 v 3 these are the internal reaction the external reaction are v 1 v 2 v 3 v 4.
So, these are the external reaction, you already know the, it is denoted by these are the external reaction and these are the internal reaction v 1 v 2 v 3 and now you construct these tachymetric matrices basically this. I assume that is metabolite A, this is metabolite B this row 1 is metabolite A, row 2 is metabolite B and row 3 is metabolite C. So, now, we can see that for metabolite A and the columns are v 1 v 2, v 3, b 1, b 2, b 3, b 4.
So, for metabolite A, we have -v 1 that is going with metabolite A is consumed. That is why it is -1 and b 1 is positive one, because it is giving metabolite A. So is b 1 is 1 and b 2 b 3 b 4 or 0 and the metabolite B if you see that v 1 is 1 and then v 2 is -1and v 3 is +1 and then b 2 is -1, because it is producing v 2. The reaction is consuming B, metabolite B that is why it is -
1 and then for C metabolite C we have 0 1 -1 0 0 -1 1 and with the reaction v 1, for metabolite C is 0 and then v 2 v 3 or 1 -1 respectively.
So, 1 for v 2 and -1 for v 3 and then we have v 1 v 2 0 v 3 v 4 are -1 and 1 respectively. So, b 3 is -1 and b 4 is +1. So, now we know how to construct the matrix and this matrix if you can enter in the MATLAB and then you define the fluxes v 1 v 2 v 3 v 4 and the under steady state dx by the metabolic constant dc/dt is 0. So, this is the steady state approximation and then you put a flux constant 0 the lower bound is 0 then upper bound is 10 for all the fluxes and you maximize the production of B said that is the flux through B 2.
So, suppose I maximize the flux through B 2. So that is the flux which is coming outside this one. So, in this problem, we want to so far, I told you, how we can actually maximize the particular reaction, for example, biomass is maximized. But here for a small network, while you do not have the biomass equation, you just define that the metabolite the flux reaction flux b 2, you want to maximize.
(Refer Slide Time: 31:29)
So, for that, you have to actually use a linear programmer, linear programming that you want to use that this is the command linprog, linprog will be used to optimize this problem. So, this is the problem definition I am just using here and then you can see that how we can actually solve. So, the I define objective function here f is your objective function A is a symmetric parameter is a dynamic symmetric when b is a metabolite concentration that is dx /dt equilibrium is the steady state parameter for S.
And then B equilibrium is basically steady state metabolite concentrations change that is dx / dt = 0 and we have the lower bound lb and upper bound ub. So, this thing you have to provide lb and ub you have to provide objective function also you have to provide basically which reaction you want to maximize or minimize that you have to specify and is basically this storage symmetric matrix idea to provide and be there is dynamic concentration that we have to provide.
So, b equilibrium is 0. So, this is 0, you already know b = 0 and A equilibrium is this tachymetric matrix. So, this thing here to provide and then you can actually maximize and minimize the objective answer.
(Refer Slide Time: 33:10)
For that we define that problem for A = S dot v = 0 and then we will just call it as S. So instead of A equilibrium we call it as S here in the previous you can saw that the A equilibrium is nothing but that is a tachymetric matrix and b equilibrium = dx / dt that is all 0. So, this one is all 0 and the lower and upper bound you can specify the lower bound as all 0 the flux vector is 0 and the upper bound is maximum value of 10 and the flux constant you have already defined likewise when there will be a problem.
So, there will be a problem says that maximize C dot v such that S dot v = 0 and v 0 to maximum say for 10 this one small v i So, this this is the maximization problem we want to define and they have defined the lower bound and upper bound. Now, we want to solve this equation in MATLAB how to do that and for that you have to define your objective function objective function again you can define by the row matrix.
So, all the elements are except b 1, so this is the position of the b 1, b 2, I think it would b 2. So, this, this one, we want to maximize this one. So, b 2, we want to maximize the production of b, that is through b 2. So that is why I have maximized b 2, that is why I kept -1 there and then why want to maximize that b 2 how much maximum b production b is produced, that I want to maximize. So, your equation for the MATLAB command is linprog f, A, b, S. So, we have to define f, you have to define A you have to define b, S dx/dt lb, ub. So, let us try.
(Refer Slide Time: 36:11)
So, this is the stoichiometric matrix, you want to use it here. So, let us define S and this is -1 0 0 1 0 0 semicolon, 1 -1 1 0 -1 0 0 is a end second row and then you enter the third row is 0 1 -1 1 0 0, -1 1. So, this is stoichiometric matrix S -1 0 0, I am just checking -1 0 0 1 0 0 0 1 -1 1 0 -1 0 0 0 1 -1 0 0 -1 1. So, this is metabolic network stoichiometric biochemical network that we entered. Now so this defined S and let us define f, your f is basically 0 0 0 -1 0 0.
So, this define your objective function, so I kept -1 here, because you want to maximize b 2 and then your dx / dt, dx / dt is basically column matrix, but all the elements are 0 this one basically, so this one I am defining, so it has 3 rows. So, it is 0 0 0 and then the lower bound and upper bound, lower bound will be equal to so we are going to use command 0s and you have 1 2 3 4 5 6 7 7 elements, we put 7 and 4 that you only put 1, so it becomes so this is your upper bound.
The upper bound have all a lower bound lower bound value are all 0, the lower value this one as well. So, I am putting a constant lower bound, all the elements are actually 0 0 0 0 0
because the lower bound is 0 here and then the upper bound will be once 1 7 and then multiplied by 10 because the maximum value is 10. So, this is become, this becomes your upper bound. So, we have a lower bound and the upper bound defined. Now we will run the command linprog.
So, this is the optimizer linear programming solver which is available in MATLAB that you can use to run the solver. So, you type this command f v fval equal to linprog f before that I define A is a empty matrix A equal to empty matrix are defined nothing is there A is empty and also b is empty so b and A are empty matrix. So, let us trying the; running the optimization.
V fval equal to linprog f A b then your tachymetric matrix S, dx/dt we are already defined, then lower bound and upper bound. So now you can see optimum solution found. So, this solver is saying that you have optimize it and found the solution found. So, the maximum value of the objective function is 10. So, you can ignore the minus sign because you choose a minus over here.
So, this says that the maximum value is 10 and the flux associated to this reaction v 1, v 2, v 3 are 0 0, so v 1 v 2 are 0 this says that your solution is 0 0 0 10. This only v 1 v 2 v 3 is carrying flux and then v 3 v 3 is 0 v 4 is 10 and b 2 is b 1 is 0, b 2 is 10 and b 3 is 0. So, you have v 3 is 10 b 2 is 10 and then b 3 is 0 and b 4 is 10. So, let us draw it what do you how much you are getting here? This is 10 that is a maximum value we have chosen and then this is 0 v 3 is 0 v 4 is 10 and v 3 is 10.
So, if the flux is going this way, then you have maximum value of v 2. So, this is one of the popper that you can this is not the linear solver you are able to actually solve the pathway that has cutting flux v 4 v 3 and v 2 that is the maximum flux we can go through v 2 then minimize also you can minimize the flux what is the minimum flux it can carry also you can calculate the null space. So null space of S this is a null space we have.
So, how do you minimize for minimization, you can actually instead of f you can put minus -f so then you can see the minimum value of b 2 is 0 the minimum value of b 2 is 0 in that case we have v 2 v 3 are 10 and v 3 v 4 are 10 whereas b 1 and b 2 are 0 and also v 1 is 0 so let us plot it in the network so we have v 1 0 and v 2 v 3 are also v 1 v 2 are also 0 so this is
when you minimize this becomes 0 and then you have v 1 is 0.
This one is 0 so that is why b 1 is also 0 v 1 is 0 and b 1 is also 0 and then v 2 is 10 and v 3 is also 10 b 4 is 10 and b 3 is also 10 is going like this now you can see that for b 2 b 0 this is the only pathway we have so this way you can minimize the objective function you can choose any function here we minimize b 2 maximize b 2 minimize or maximize b 2 but you can choose any one of them and in flux you can consider for maximization or minimization.
And you can run it in MATLAB and know which pathway it is following. So this is a simple network optimization we have done in MATLAB using linprog as the solver which is freely available in MATLAB. So hopefully we will be able to run MATLAB easily on your computer unable to calculate this value using the network. Any network you can optimize and find out what is the maximum flux going through the reaction what is the minimum flow going through the reaction you can easily calculate.
(Refer Slide Time: 48:14)
So we are next the next thing we learn about how to optimize the E.coli core model where we have many reaction that we want to optimize so in the E.coli reaction we will consider only the coal metabolic network that is different amino acid formation and then glycolysis dcs cycle and those networks you can optimize and find out what is the maximum biomass it is flowing how much is the growth rate all those things you can calculate so in the next class we will do the E.coli core metabolic modeling optimization.
(Refer Slide Time: 48:59)
(Refer Slide Time: 49:00)
Then these are the references you can follow thank you. Thank you for listening.
Log in to save your progress and obtain a certificate in Alison’s free Diploma in Metabolic Engineering online course
Sign up to save your progress and obtain a certificate in Alison’s free Diploma in Metabolic Engineering online course
Please enter you email address and we will mail you a link to reset your password.