 I'm very happy and glad that I can be here after seven years. So I actually, an Anumi of this spring college, like I attended the 2016 versions. And after seven years, I have time to come back here and just do a little bit of tutorial on what I am interested in now. And so, as you have seen over the last two days in CUNY lectures, so we want to make a physics-style model for some biological systems. And so, CUNY will give you a broad overview of his model. He has been developing over many years. Today, I will just focus on one of his model, which is a bit different from the model you have seen during the last two lectures. And so, let's go to the model. So this is actually the paper he published, like in 2007, on the evolution of a gene regulatory network in which he considered two effects. So one is mutations, no effects. So mutate by mutational effect, which means that you have the chain in the genomes and another effect from the noise in the environment. So basically, the main questions, the main question is that if you have noise which represent environmental effect, and if you have the genetics, let's say it's a mutational effect, so this is environmental effect, if you have these two components and both of them can act simultaneously on a biological system, then depending on the interplay between the first and the two, the system will display some emergent behavior. And that emergent behavior is something that you can interpret as a phenotype. So phenotype in this particular context of that model, for instance, he understood as the gene expression level. So you have a bunch of gene that are connected in some network and then for each of the gene, you have a state variable Xi for simplicity, you can make it between zero and one. And then the state of this gene regulatory network is characterized by this vector of n component when n is a number of genes. So this vector, the configuration of the genes is a phenotype in this context and the specific network structure and so not only the network structure, but the strength of the coupling constant between the genes understood as a genotype. So the basic idea is the following. So if you have these two components system, so the first is genotype and the second phenotype, how do they affect each other? And in the more traditional setting, what is considered is that you fix the genotype which means that you fix the network configuration, you don't change it and you only observe the dynamic of the phenotype. So you will have some dynamical system corresponding to the dynamic of the phenotype. And so the new component coming into play from Kuni model is he considered that phenotype also involved, but on a slow time scale. So this dynamic is on the fast time scale. So you have a couple of the two dynamics happen on different time scale and one will affect the other. So you create a feedback loop between the two components. You have the genotype affect the phenotype dynamic and phenotype on its own will affect on so the genotype. So you have a feedback loop between the two which make the so called notation in the literature called evolutions. So the system that have that type structure like you have one type degree of freedom and you have another type degree of freedom. And if they have the dynamic couple in such a way that the evolution of one affect the other then it will be called co-evolutionary system. And in that type of system which you go beyond the standard dynamical system framework you can observe a lot of interesting behavior. So say again, my task is not to try to do some tutorial to show you how to make a simulation of this type model but more or less like inspired by Kuni style. We will try to go through some basic concept and how to formalize this structure, right? So how to formalize this type co-evolution system. And I already send a code to the secretary of the workshop of the spring college. So after ending this section you can just get the code and then work through it yourself. Of course I've been so a bit how the code structure and the basic behavior during the tutorial section but the first idea to get the general picture of what we are trying to do. So okay, now come back to the first step. Okay, yeah? It's there. The sources of the noise is some way for simplicity as in physics. You can consider it as a term of fluctuation. So you couple the system to some heat bath then it produce some fluctuation. So this is the noise. So this is why I call it environmental effect. But so this is your system and this is environment. The system is coupled to the environment but it also have some ancient sick property. So the ancient sick property in this context is a genetic network. So the network structures are coupled for different genes. So the question to be asked and to be answered here is what is the interplay between the external condition and the ancient sick internal property of a system and how it's been resolved in the emergent behavior in the center of physics. So connective motion like in condensed matter physics. Yeah? Then the other question. So this is a general picture before we go to the specific detail how to model this system, maybe how to run the code of blah, blah. But this is a thing that I am trying to show you. So now let's go to this specific part. So you have two components, right? You remember that you have the phenotype and you have the genotype. So first we consider the phenotype. So the phenotype is presented by this vector of n component, right? So the one that I brought over there, you have a vector of n component, which means that you want to have a map. That shows you just how the system is born over time, right? And so just today, more or less, we need to show the deterministic system. But in this model, you will consider some effect of noise as you see what we are interested in. And starting from this general structure, one specification, maybe I wrote it out here so that you can see it better. This type equation, you will see it quite a lot in many different contexts. I first write down, then I say what, in which context and in which type model you can observe it. So you have the derivative, the time derivative of the gene I, right? Gene I just want components, yeah? And here, you will have the decay. Every one of you have seen it in the previous lecture, right? Why do you need that? Just to make the system about it. If you don't have that decay, everything go to unbounded regime. And then you have this tangent hyperbolic function of something. And then you have the noise. So first of all, why is this tangent hyperbolic? Actually, any sigmoid function should work. So by sigmoid function, I mean the function that go from zero to one here. Even x can go to infinity, right? And the tangent hyperbolic function is just something go from minus one to one. So, as a physicist, you love the symmetry. So this is why you want to have minus one and one, for instance. So this is just for simplicity. Actually, the hill function today you see can also be related to that one. So anyway, this is just to give you some argument why the dynamic is started by this function. Now we need to specify the argument, right? Any of you have some idea of what is the argument inside that should be? Remember the question you want to see. How the effect of noise and system property affects the dynamic of the genes. So here you see, I showed you the noise. Here is the decay. So do you may guess? This need to be the system internal property, right? So, and now since you specify that the system, the component depend on each other by this network. So this need to be a function of the graph representing the interaction between the genes, right? This is, would you agree with that? Argument, the argument here, yeah? And the simplest way as in physics, you try to capture the common reasoning is that only pair-wide interaction matter. Of course, you can consider high order interaction, the interaction between many components, simplex or complexes, whatever, hyper edges. But for simplicity, you can just consider the pair-wide interaction between two genes. Then here, the sum, you can represent that kind of pair-wide interaction by the sum over on the node j in the neighborhood of the node i. So this is my notation. So delta i represent the set of neighbor of i and j belong to that set of neighbor. And here is the way you represent the effect of the network structure on the state of the node i. So you see that this tangent is just what people know in the context of neural network activation function, transfer function. So for those who know that already, so this will be the dynamical equation for one single node. Now, for making it a bit easier and chuckable, you can specify the noise here with some y noise, right? So you can specify some psi i and that's the y noise, psi y, psi t prize and here two sigma square and delta t minus t prize. Something like that. So you consider the y noise for simplicity, but this will be the dynamical equations that you are interested in to study the dynamical evolution of this phenotype. So for given network, the network just specifies the inhibitory and promoting interaction between the gene. What happened with the state of the gene? So this is basically the first step that we want to model, the dynamic of the first component of the system, the phenotype. You have any questions? Do you know it here, not velocity? It is a step here. It is. Ah, sorry, sorry, sorry. Okay, okay. I might still confuse. Sorry about that. Sometimes I denote with a dot and sometimes I, yeah. Sorry. Yeah. Yeah. Thank you for correcting me, yeah? Hi, I never have a proper discipline of networks and so could you just talk a little bit more of the hyperbolic tangent and what is inside? You said it comes from networks and so on. Oh, okay. I'm sorry. If I don't make it clear enough, you have two factors that can affect the dynamic of the phenotype. So this is a noise, you see it. And the other need to come from the interaction between all the gene in the network, right? But of course in the network, these particular genes are not connected to all other genes. It's only connected to a subset of the genes in the network. So this is a notation for the subset denoting the neighbor of i and then all the neighbor of i can affect i, for instance. And this is how you represent the effect of the genetic component of the network structure on the dynamic, okay? So far so good. Okay, so if there is no question, now we need to a little bit more specific. So in the disorder system, so the topic of my most interest, people normally consider different ensemble for that main matrix. So it can be a Gaussian ensemble, for instance. So every J ij is distributed according to some Gaussian distribution with zero mean and some variances. So as a starting point, you can try to simulate the dynamic of this system with the assumption that the network effects, so the topology of the network effects and the value of the coupling strength is distributed according to a Gaussian distribution, right? This is first step, now. Okay, now I can come back to this picture, okay? It does not work anymore, wow. Okay, now it works, sorry. Oh, yeah? So why did you took a term like that? Why is the argument is like that, like J ij x i, why it's not something other, x i to the power alpha or something else? Sorry, sorry? Why did you took this term like that, argument? Like, would it be like something x i to the power alpha? Higher order terms in x i, would it be something like that? Why we don't take the higher order, x cubic, square something, right? Yeah, yeah. That's very simple. Remember what I said, it's zero to one. Or maybe in this case, better to see it is minus one to one. Then if you take the higher power, right? It just becomes smaller and smaller. And by the way, it is very similar to the structure that you know in statistical physics. Like in easy model you see for instance. So you consider the linear term or cost because the non-linear term we had to do it. And now I thank to you, I forgot one thing that I want to mention. You can put a factor beta here, similar to the role of the inverse temperature in statistical physics. So this will be the final equation with the beta here. You of course can set it to some number. It just rescaling accordingly to the strength of this term of the Jij. But if you want to make it explicit, you can write explicit here. What would be beta in this case? So you can just write beta here as a parameter. It's also five, yeah? Thank you. So these x's represent your phenotypic values, right? Yeah, the x represents the gene expression level. And you have many genes. So each of the gene has its value, xi here. But to specify the full phenotype, you need to specify the n dimensional vector x. Okay, so this equation you need to understand it's the equation for instance. If I have 100 gene, you will have 100 of such equations. Okay, so this. And if you have 1000 gene, of course it's to go to 1000. So this vector is a representation of the phenotype and each element of this vector is one gene, the expression of one gene. Yeah, the expression of one gene, yeah. Okay. And you were talking about this interaction between genotype and phenotype or the co-evolution picture that you were talking about. This is a general picture. We have not reached that point. Ah, okay. We now specify again fixed genotype. So this is the first step. We don't consider the co-evolution yet. Okay. Yeah? Yeah, thank you. So just to check if I understand, if there were no relations between the genes, you would have like an additive gene thing, like an additive gene model. No, if this is zero. Yes, if this is zero, what happens? If it's just zero, then you just have minus second noise. Okay, and what would that mean? This is a standard long-term equation. You have the diffusion, Brownian motion pattern, the return, the return to what zero? But what would it mean biologically? I understand that, but biologically, what would that mean if you don't have the gene relations? Okay. The genes would just- This just means simply something like that. If there is no interaction, there's no way to make the system order. So you know the central thing in physics, where you want to describe the transition between order and disorder phase. No interaction simply means just decay to the zero state. Everything just becomes, so let's put it this way. So if on the x go to zero, this is nothing that happened interestingly, right? Because you want to see whether the gene can be on and off, and on is plus one and off is minus one. So if everything is zero, it's just how to call it, the sleeping state, or anyone have the better name for that? If on the x is just zero, that state or sleeping state, dormant state? Yes, okay. What you would like to call, just like a paramagnetic phase. So it's not an interesting behavior. So this is why you want to see the effect of interaction. I'm sorry, if I am not good biologist to explain your question. So does it make sense you are trying to find the transition between phenotype to genotype? So is it what you are going to do? You are trying to find the transition between phenotype to genotype? No, no, no, no. Wait a second, we go to that slide. First we need to understand at least the modeling framework. Okay. So it's just the time evolution of the phenotype we are now studying? Did I answer your question or still? Okay, okay. So can I continue now? Can I move on? Okay, so we will do just a simple example. If I stimulate that dynamic, right? What kind of behavior can I observe? Let me take a few examples for you. This is with the evolution. So as you guess, okay, what happened here? Okay, so this may be interesting for you because it's related to George's questions. So on the left-hand side are three representative trajectory of three random genes, so one, five, and 11. And on the right-hand side is the histogram, the distribution of the value of the genes at the final time. So this is the left and this is the right-hand side picture. And for a particular choice of beta is more enough. And sigma noise strength here is also kind of small. And you see that thing just decay to zero. Think decay to zero, the state may be just great to do a question. But when you increase a bit the noise one instant way, if you increase a bit the noise, still things is not much interesting. It's still fluctuate around zero. It just have higher variances. And now if you consider the strong enough beta, so the inverse temperature, then you can have this guy picture for instance. Now you see it's the more interesting way. You have some component go to one, so it's one here go to one, which means that it's steady on, this genie on, starting after some transient time. But for other genes, it can behave like limit cycle. So it just have not settled down in some units fixed by a tractor. And of course, you can see nice there with you the code. You can play around with different set of parameters and see more complicated behavior. But basically there will be three regimes. One is a zero state, everything just decay to zero, right? At the small beta as you seen in the first picture here, where is that first picture? Now you will just see everything just decay to zero. And you will see some more interesting state which genes can turn down. And some of them still have persistent like limit cycle like behavior. And as a state in which only just go either to on and off. So, and so this is the basic, okay, yeah? My question has to do about the fitness function because I'm trying to understand the... Now I will explain to that. First, we just think it's simply. Yes, I don't make sure to do. When I write a fitness here, basically it's a little bit of a, I abuse a bit of the notations. I should call it just the final state, the state of the genes at the final time point. So you just basically plot. So you have this vector, right? X at the final, in that picture it just at the time 400. Then you plot the distribution of the gene at that time. Of course, you can run it longer. Every time you run it longer and you observe a different time window, you will see slightly changes. But this picture will be something established on the system and to the stationary state. Yeah? So far so good, no? More questions? So how I pick the sigma? This one, this is the component in the equation. This is the sigma. This is the noise, it specifies the correlation of the y-noy term. So for this particular equation, you have two parameters. One is the inverse temperature and the other is the noise term. Of course, I don't want to go into detail. When you write a code for that system, you have to specify a bit more. So you need to, for instance, specify the link density. So you want to make the network dense or spot, and then you generate the network of gene regulation. As I said, it's typical to be chose as a random Gaussian distribution. So once you have this network specified and you specify the beta and sigma, then basically you just write two lines of code to simulate the ODE. So it's not that complicated, right? But the behavior of this super simple model, so the model which basically I've written, have written it here, can generate a very different dynamical behavior. So just a command. Again, if this j's, it's listed as a branch random variable. So for those who know the terminologies, quench means the coupling constant here affects over time. It does not change over the course of the evolution of the genes. So if the coupling is quench, then there is a theoretical framework, the so-called dynamical mean fin theory, which has been developed to make some analytical statement regarding to this equation. So for those who want to know it, you can check some old paper by Sumpolinsky or more recent review by John Hed or Tom Cullen. So if you are interested in solving that equation, just come to me, I can try to use some reference. But this, now it's just the first part of the model. You remember that we have the second component, right? The genotype. And this component will be very important because the behavior I showed you here for given network of genes interaction will be changed quite dramatically when you have the evolution of the genotype. By the way, it's a common for the mathematical treatment to assume that the genotype involved on slower time scale so that you can do some kind of asia-abatic enumeration or how they call it, pastel allulink, right? How they call it in the spin glass literature, so called pastel allulink approach. But in this particular model, CUNY also assumed that the evolution of the genotype happened on the slow time scale. But now this is not for the sake of mathematical treatment anymore, but it's really reflect the reality because in reality, you expect that the genome take million of years to inform. So it happened much slow on a much slower time scale than the genotype. And now we go to the second step, how to model the evolution of the genotype. Is there any question? Because this second part will be more interesting for you and I just want to make sure that everyone of you on the same page, yeah? Okay, now let's go to the second part. Which, okay, now let's go to the second part. So now I will try to write it in the algorithmic manner so that after reading through, you can come up with the old code, no? So, let's call it step one. In the step one, you simulate this system of ODE, right? To get the state of the phenotype. Now at the end of this step one, we go to the step two, and now to answer your question, we need to define a fitness. Because what is a fitness? So fitness in biological understanding, it's just proportional to the reproduction rate, but you can think about it as a function of the phenotype. So if you interpret the fitness as a function of the phenotype, so the fitness need to depend on this vector x, right? Yeah, because fitness is based on the observable chat of the system and then it's a function of this observable chat. This is x and now everything boils out to how you specify this function, right? If you make that function too complicated, then you may get into trouble. So you want to make that function simple enough. So again, now in physics, if you have a state vector, what is the simplest function that you can have? Which function? It just, maybe you want to make it like that, right? Like the marketization, right? You just sum over all the components, right? And then you divide by that component, the number of components. So basically, you consider the mean activities of the genes as a fitness function, right? And there is a little more than that in Kuni model. It's the sense that actually he divide the set of the gene into target and non-target. And non-target. So to make things maybe more interesting because his assumption is that biological system consists of two parts. One is a functional part, so which is responsible for performing some function and some, you can call it HEP assistant part, which does not play any role in determining the fitness. So maybe here, the sum will not own the node but over the subset of the node belongs to this target set of the gene, yeah? But anyway, it does not matter. It's just a very specific way to model it. You can just think about it as a mean marketization for given set of genes. So now this is a step to define the fitness based on the steady state behavior of the gene dynamic. So here I forgot to say that mathematically speaking, you want to consider the X-infinite. So X-infinite is a limit of T go to infinity of that X, which involved according to that stochastic dynamic. But in simulation, of course, you cannot run the simulation infinitely long. So actually this is just a T after some transient time. So you run the dynamic long enough. You assume that the system settles down into the steady state and then you start measuring the phenotype and then you compute the fitness. Are you okay with that? No questions, right? Yeah. Yeah? That's for each gene, I'm gene number 8H5. Yeah? And the eye is all over time? No, no, no, no, no. Eye here is the index, the label of the gene. I am gene number one, you're number eight, yeah? We sum each other. Yeah, we sum each other. We form a collective. So we want to think about emergent behavior, right? Collective behavior of a large set of degree of freedom. The gene state that you got from here. Hi, how can we define the target cells inside Big D? So it's very interesting question. It's just very important. I also asked Kuni many times, how do we distinguish target from non-target? Actually, in this type abstract framework, you can assign so gene number one target, gene number two and non-target. So this is a choice. But maybe in the realistic setting, you need to have a way to distinguish between target and non-target. So so far in his model, he just predefined certain number of gene and their indexes. They are the target genes. The rest, the remaining is non-target, yeah? Okay, I can continue, right? So we specify the fitness. Now, as in biologies, you have the selection process. So selection process is a process acting on choosing the best fit in the population, right? So how did he model this selection process? Now it's been more tricky. Where can I clean it? Okay, let me clean it here. Ta-da, ta-da. So let's put it this way. At T9-0, you have a, let's put it this way, you have a five different network. So you have five different genotypes, right? Each of that have one network, yeah? Each of the individual in the population have one network, that is the raw as a random instant from the ensemble of random matrices that specify by the random Gaussian distribution, right? So you have five network at the initial time and then you run the dynamic. And then at the end of the day, you compute the fitness for each of that network, yeah? So it gives you F1, F2, F3, F4 and F5, yeah? You have four different fitness functions corresponding to five different fitness functions corresponding to, sorry, five different value corresponding to five different network. And then you just rank them, yeah? You need to rank them according to some order. So for instance, if you observe F1 smaller than F2, and maybe smaller than F5, then F4, and then F3 is the maximum. Now for simplicity, GUNI specifies something like only the first, like the two most highest fitness network can survive, can reproduce to the next time step. So this is how this is chosen among all the network. So you choose two among five, which means that the selection pressure is like 40%, right? You choose two network among these five network and this is the step three. You choose NS network among the total M network. So M equals five here and NS equals two here. You choose two among the five network and then you do random notation. So the random mutation is the following simple recipe. So maybe many of you have no system keyword again beforehand, rewiring. So rewiring just means that this is GUNI and this is GUNJ. So mutation happens in a way that you remove the connection between INJ and you flip it to J2K, for instance. You do it randomly. So you choose a pair of GIN, random, and then you create a link that had not existed before also at random. So this you perform on the last two network here, on the two networks that have the highest fitness value. This clear? Because this is a crucial step, otherwise the dynamic is just simply simulating this evolution of phenotype without alternating. The network structure in the population. Yeah, do you do the mutation step only for the top fitness, top fitness network, yeah? No, no, no, don't worry. So what you got till here was like some kind of optimization of the some things and then once you get like let's say the best to colonists, now you mutate like a higher level thing. Yeah, you can think about that way. So you have a system which operates on two level. The low level is a level of phenotype and the higher level is a level of genotype. And on the level of genotype, what you operate as what you said, you try to maximize his fitness function. So this is how you retain only the top fitness network among all the network that you simulated at the beginning. Thank you. Yeah, and once you make this random mutation, sorry. So I'm just gonna clarify about the rewiring. So it seems like there is a lot of parameters that would be involved, say how many connections would you rewire or like do you want to define those? So one thing again just because the assumption the underlying assumption is evolution happened on a very slow time scale, which also means that you typically assume the mutation rate is very small. So to draw question, which means that at any given time step, any given generation of the network, what you want to mutate, you rewire only one link. So you kill one link and you create one link. And this is how you make the new network into the next generation. Okay, so I'm a little bit confused about what the rewiring means in terms of the changes in like my values, for example. So I have a value for each node, each node has a value. I rewire and what does it mean for them? You don't rewire the value of the node. You rewire the connection between them. Okay, so I rewire like the term that has the... So it's the same like in the social media. Instead of following chump, you follow me. And for the next generation, I have the rewire. Yeah, I think... Just to recap, you have the equal number of networks as the number of phenotypes, is that correct? Yeah. And then for each one of them has their own dynamic. Yeah, each of them you have their own dynamic and at the end you compute the fitness. Okay, and for the top fitness, you take those network and rewire them a bit and take it to the next generation. Yeah. Okay, thank you. Any other questions? And just to be sure, this rewiring process is random? This is random. You pick a pair of nodes at random. And you remove the interaction between them and then you randomly choose another pair of nodes and then you create link between them. Okay, and how many connections do we change for how long generation? One. Only one. Only one, most of the time only one. Okay, no more questions? Okay, so basically this go through three steps, right? One, you finish a step, three, you go back one. You make that cycle many times. So each generation consists of three steps. Just repeat it and you see what comes out from that dynamic way. Now let's see what can we find. Okay, let's go with the pollution. Sorry that I forgot I need also to specify two additional that go into observable to make the thing more interesting. So, where can I go? Okay, let's do it here. So you have this phenotype, right? It's important according to this stochastic dynamic, yeah? So even for a given network, for a given network, if you run this dynamics, you simulate it every time you have slightly different random trajectory, right? So for a given network, actually you can define the mean value of the X as a function of the network. So let's put it this way, let's put it with a bar, maybe better notation, and then to specify some distribution which is well to the X and then of the value of the fitness of the X, for instance, you will take it as a function of the G. So this one is a function of the network. So basically you consider the mean of the X for given network, G here stands for network, right? Yeah, it's just a quantity. And you can do it for given node i. So you can marginalize it, yeah? So the distribution here is over what? In the sense that the idea is that you have multiple the distribution here is generated by this random ensemble of trajectory. Looking at ideally for a very long time. So the idea is that you have alternative stable state. You can have multiple attractors as I showed you in the first place. So you have a distribution of the state of the node i in the stationary phases because this is inherent disorder here. You can end up having multiple attractors cases. So under the assumption that the system does not go to chaotic regions and in the present multiple attractors, then you can specify that distribution for the state of a G and i for given network. So this is one quantity. And then you can consider the variant of that quantity. So the variant of that quantity X i, for instance, yeah? The variant of this you denote it as a VIP. So VIP here is what the Coonico needs. Isogenic variance. So for given network, you run the stochastic dynamic and you end up having in different state. And then you take the ensemble of reach over the set of state. And then you can consider the variance, right? In that ensemble. And this is the first sort of randomness. Now you look at this recipe. You see that you have the updating rule because of the rewiring for the genes. So which means you have second sort of randomness. Naturally, for a second sort of randomness, you would expect to see another type variances. And this is what he calls mutational variance. So the mutational variance is first need to be defined as a following. So once you have this, the iconic IG, right? Here, XIG with that. And then you consider the distribution of the gene. So the distribution over the network ensemble, right? And then to specify a measure in that space, this will give you another vector, another value, you call it XI. Now this is the average over the average. So this is the average you get for given network. And then, because you have an ensemble of network, you can define a second average, which is this guy. And of course, corresponding to this guy, you will have a variance. And this variance is the one tactical mutational variance. I don't know what level of detail that you expect. If you want, you can talk more in detail about that. Otherwise, I will just go quite fast into the main conclusion of the paper and how we get to that main conclusion. Are you clear how do I define this fit of the first function? Okay, just, yeah, yeah. Of course, in order to compute that, in order to compute that, you take into consideration the time scale separation or something. So once you set that the time scale is different between the dynamic of the phenotype and genotype, they effectively, during the cost of evolution of the phenotype, the genotype are fixed. Since they are fixed, there is only one randomness now, which is the randomness play in this equation. So that is how you compute the first quantity, yeah? And now, since you have the evolution, which is basically the mutation of the network, yeah? Then even you start with the same phenotype after you get into this second step. Because of the law here, right? The random rewiring law, you will get to different state and then you need to compute the variance of the steady state when you perform that random rewiring. So this is the second quantity here, yeah? This is confusing? I'm sorry about that. This is the first time I explain that this may be not clear for you because of the way I explained it. So I'm happy to repeat and explain any questions. I can go through the code, but now I just want to emphasize the physics and the main question first. Because otherwise, as I said, you get a copy of it through the email and then you can play at the old page, and it's better, right? You can look through every single line, figure out what is reasoning behind that. But the main idea is, again, we want to consider these two effects, the noise effect and the mutational effect on the dynamic of this gene regulatory network. And okay, so you have two variances. This is a VIP, it's a variant of that guy, and okay, I need to clean something. So what is that one I need to clean? And you have the variant of that HIG here. So that quantity is what you can call it VG. So you denote it VG, and here you denote it VIP. And here is a variant of XI for a given gene. Okay, so you have two variances, yeah? And the most important conclusion in CUNY paper is that the relationship exists between the two. Actually, this is an inequality. So this can go both sides. So do you expect that there is some free transition when the two equal, the two quantity can become equal at some critical value of the noise, sigma, but for small noise, you will have VIP smaller than VG. And when you increase the noise, two things become equal. And if you keep increasing the noise, VIP become larger than VG. So this is the main conclusion of that model. Of course, I am sorry that I'm not good at biology to explain to you all the, how to say, biological relevant of that inequality, but to, shortly speaking, we just mean that you have a constraint on the evolution of the system because it's a relation between two level dynamics. One operated the phenotype level and another on the genotype level. So depending on the strength of the noise, you have some different relationship between the two dynamics. And it's reflected in this inequality. Now I will show some behavior and then we go to the code for the whole night, yeah? Okay, so let's stick with this first example. So small noise, yeah? As I told you, small noise, you will have these mutational variances, VG, larger than VIP. So this is what observed by simulating this dynamic. Yeah? And if you increase the bit of noise, let increases, what did you observe? So this is noise equal 0.1. If you increase this to 0.5, you still see the same inequality whole. So the VIP, the isogenic variant, is smaller than VG, the mutational variants. But then, once you upload the critical value here, seem to be sigma equal to one in this state, right? You see that the two variants become more or less statistically equal to each other, yeah? And if you keep increasing the noise, now you observe that this whole, so VIP will become larger than VG's beyond some critical level of the noise. So just a demonstration of this inequality, how it's work in practice, yeah, question? Please, could you give us some intuition about why the VIP has a linear relations with the noise? So intuitive explanation of that? Like if I get you, you mean the VIP increase as we increase the noise? Yeah, this is true. So the VIP, of course, according to that dynamic, sorry, I just maybe rewrite the curve. So yeah, this is true. This is a tangent hyperbolic here and then the noise term, which proposed not to sigma. So the VIP as you observe, it should always increase with increasing the noise strength. It's true. So my question is, can you give us some intuitions why this is possible? So the one way to think about it, let's think about it in the physical picture. So it's not the correct metaphor to map the system to a spin-glass system because in the spin-glass system, you have the symmetric interaction, Jij, and then you can define the Leoponopax polon and then you can like map that picture to some evolution in some free energy landscape and then you can explain it quite intuitively, like you can end up in multiple ejector corresponding to have the phase space broken into multiple erotic component. But here, the Jij is asymmetrics. So it does not work the picture of the Leoponopax polon. So at least you can think of and at least if you still remember what I showed you in the first part of the tutorial with that for the small noise case, the system can end up in multiple ejectors and each of them can be think of as some kind of metastable long-lived state and if the system just with the initial condition turns out to be in the basin of ejection of that metastable state, it gets into that. And because of that, if you consider this VG, then there is a huge differences between different ejector. So if the system in one ejector and you compare with the system in another ejector, you can observe the tremendous variances between different ejector here, even for the small noise case. This is what really you observe when you simulate this dynamic. But then when you increase the noise, actually the noise act in the beneficial way. So normally noise is associated to some harmful effect to the system. But here the noise on the opposite, it can be kind of resemble to the so-called stochastic resonance system in the sense that the noise actually stabilize the system and help you to, so if you think of this in terms of some landscape picture, the noise just sharpen and make the landscape become more full now. So in the case with small noise, there are multiple local minima in the landscape picture. But in the case with noise increase beyond some critical value, one among those peaks get deepens and then this will be the most probable state of the system and it's redo the variances across a network, yeah? No, I don't understand what you said about the multiple attractors because if you go to the small noise case, the isogenic variances very small. No, no, no, no. Actually, correct me if I am wrong because you are more expert on ecosystem. So in a bunch of paper recently by Tobias Gala or Bulleen or whoever, Pankai, they saw the dynamic with all the noise of the low-carbon model, generally low-carbon model, for instance, for specification. Do you actually observe multiple attractors? That I understand, but what I mean is that if I look at the picture in the case of Sigma very small, you see that the isogenic variance is not these pictures, the one for Sigma equal to 0.01. Let me go back to that, 0.01, where is that, 0.01? Yeah, or this one, the variance is extremely small. So it looks like you don't have multiple attractors in this case, right? No, it's really small, right? So you have one attractors and then there are perturbations around these attractors. Yeah, this is the more random. Exactly. No, mutation, with respect to that attractor, but actually it can also be so that in some given value of beta, it can end up into multiple attractor and limit cycle as I saw in the beginning, yeah? So the idea is that for very small noise, the system can end up in very different type of dynamic. And then this is the observation that the VG is larger than VIP. It's quite counterintuitive. When I first talked to Gudi, I also don't make sense of that. How can it be a small noise deterministic system, one single attractor? Why in small noise case you have multiple attractor? I don't get that. So it's also my technique quite a while to understand that, to appreciate that fact, you can simulate that dynamic and you can figure it out. And it's also written in the paper that you want to have some, like in the picture of what you just asked. So the dynamic in the small noise and the large noise can be mapped into different type landscape picture, even though, say again, there's no way to define the function product system with asymmetric coupling. So, yeah? And the other question is the ensemble of the networks is basically the distribution of the networks within a population which is evolving or is like a random drawing of the networks? Only at T-zero. At T-zero, which is Gaussian random ensemble. Okay. But then, after the T-zero, every time, at every generation, you mutate the system according to that fitness selection and then it will become not random anymore. No, no, yeah, that's my question. So when you calculate the VG in the picture, you are looking at the VG of an evolved population? Of the inborn population, yeah, of the inborn population. So this is why, actually, you see in that picture, the VGs, it plotted also function of the generation. So if I just specify the network at T-zero and I fix that, it makes no sense to plot this kind of picture. So here should be right as a function of time to make it clearer. I guess just a bit of clarification. Why are we, since, for your example, you have five networks and you chose a top two? Wouldn't say the number of two networks matter like in the long run? And at the same time, why the first two, not the last two? Why are we choosing like the first? We choose the first two because it's one which has the highest fitness. So the biological evolution happened in the way that only the fittest species survive, the one which has lower fitness extinct. I see, so yeah, so most of the time, say for example, for genetic algorithms, those that are in the last parts are the ones that are mutated with their genes. So like it's just new to me that the first parts are mutated. So back to the first question. How do we select two? How do we select the number of networks that need to be mutated? So actually, the model behind we acquire robots with respect to changing of that number. So you can choose two, you can choose four for instance. Of course, if you choose many, if you keep many network with low fitness value, then you will have to brought the distribution of the fitness. But qualitatively, the model is so that it's this inequality whole and it's not that too sensitive with respect to the exact number that is specified here. Is this true, right, Kuni? This is two weeks maybe. So if you select only nine from 10, maybe the evolution may be slow or it may not work so well. And so it depends on the selection pressure. But as long as the selection works, maybe this behavior. So I have a question. Do these results is still true if we change the fitness function? Cause I guess that even in biology, fitness can be something kind of tricky to define. So if we just invent a different function for our fitness, for the way we are measuring it, does we get to the same result or can be a completely different thing? Oh, thank you. It's very good questions. Actually, I asked the same question to Kuni some time ago because I also, I love simplicity. So of course, the marketization is a simplest way to compute the fitness. It actually turns out, according to my recent calculation that any monotonically increasing function of the x will work the same way. So no matter you have the sum over the x i and you have the half of x i, which half of x i is some monotonically increasing function, you will have the same behavior. Of course, if you take the parabola function or some weird function, the behavior may disappear. It's a little counterintuitive to me that you are making the fittest phenotypes mutate because if a phenotype or, yeah, okay, if a phenotype is, or if a genotype is the fittest, it means, at least intuitively, it means that all the other genotypes on the fitness landscape should tend towards this fittest genotype value, right? But what you are doing in this model is you're taking this fittest genotype and you're making that genotype mutate to a different value. So that is not making a lot of sense to me. Could you please explain that? Oh, this is a good question. I'm afraid that I cannot answer it correctly. So the point is that maybe it's important to keep updating yourself, right? So that your population always gets better and better. Like, I think this quickly remembers some paper, like, human being gets smarter every 20 years. So even the fittest can go through the selection and improve it, no? Okay, could you even answer your question? Yeah, of course. So even after, maybe we are a little bit, but still mutation occurs for the next generation. So the mutation still remains so there is some genetic change. And, but when VG is small, that means even if you change a little bit, genetic mutate by genes, by mutation, if VG is very small, then the change is very small. So even if you have a fittest one, then if you put some mutation, it's still basically fitting. So that means so if the system goes to small VG, then that is very robust to have mutation change. So that is good for bio-systems. Well, you were asking, no? The thing they are doing when they are killing all the other ones, they are saying like, that is the success for the population, just takes over. They are like skipping a step to do it more fast. Yeah, is there any other questions? Okay, from now, right? Now I will spend the last 15 minutes just to show you the code. I think Kuni will explain this model, I think at some other point in his eight or nine lecture, and then you can ask him more specific biologically related question. It's not my task, sorry. Let's see where is the code that I start. Okay, so this is the code. So first I told you at the beginning of the session, right? You need to have a simple ODE integration for the dynamic of the gene, right? You have the decay term, the noise, and the tangent hyperbolic function. And then you have some other function to capture the step two and step three of dynamism. Sorry, I just have a question. Like what is the physical interpretation of the transition VG, VIP, I don't understand that. The meaning? Yeah. So the meaning is it captures two different dynamical behavior. In one case, when the VIP is smaller than VG here, you can have multiple eye chapter and the system can get chopped into low fitness state. And when you have the opposite way, VIP larger than VG, actually the system get to much higher fitness state. And the system actually in the more detailed discussion in CUNY paper, this is the way he explains the emergence of robinus. Let me go to that, where is the paper again? Ta-da. Sorry, okay, is this this one? Yeah. So at the beginning I mentioned that what we are interested in is the relationship between the noise effect and the genetic effect, right? So the VIP represent the effect of the noise and the VG represent the effect of the genetic mutations. So you want to understand the relationship between these two effects. And this is the conclusion that the raw from the model that depending on the strength of the noise, they can have different relationships. And these different relationship emerge under different dynamical behaviors. Yeah? This is okay? Okay, so now it's back to this part. So the first step is to simulate this equation for the gene expression dynamic, right? And then you need to have the channel of code to simulate this step two and step three. So basically as I said at the beginning, you need to create a random assemble of Jij magic sis. And then you just do a bit of how to say specification according to that picture, you specify the number of the network that you want to simulate. And then for each of them, you go back to the step one. You simulate the phenotype dynamic and then you go to the next step of computing the fitness. And then you choose the best fit. So for each of the network among the top fitness in the population, you do the mutation. So you copy all the connection from their parents and then you rely with some probability to make the new network. And this done is no magic here. You can just try to do it. I stare at the code and then you play that. If you have some questions, come back to me. If you want to know more about maybe analytical argument or some literature on that, I'm happy to answer. By the way, to your question, Diacopo, I think that the picture makes a wrong impression here because even in this case, this is too small. Okay, now let's go back to that P equal one here. Actually, in the orange line here, you divide by factor one over N because you don't compute exactly XI but you will compute XI over N. So this factor subred the amplitude engine. Do you think that is almost no fluctuation? Actually, there is fluctuation there between different value of the XI. So just to answer your question, that maybe it's not clear why it seems like almost zero fluctuation. There is a fluctuation there. What I meant is that if you have multiple attractors, you should expect to have a finite VIP as you go to sigma, as sigma goes to zero plus. Yes, this is true but as I want to say, it's because of the, how to say, because you try to represent two quantity on the same scale and which have different amplitude and then it just makes the impression that it's zero but it's actually not zero. It's a, you can go back to the, if you don't have the evolution, so you still can see the fluctuation there. If you don't have the evolution, so this is no evolution, okay? And for the same value of the sigma zero to one, you still have the fluctuation, of course. You have the random trajectory induced by this noise. So okay, so I think I can finish, no? Any other question, no? Okay, thanks a lot.