Archive for the ‘dynamical systems’ Category

Second Bonn conference on mathematical life sciences

March 20, 2026

I have just attended the conference mentioned in the title of this post. My general impression is that mathematics as I understand it is being more and more excluded by on the one hand masses of data and on the other hand by reliance on computers related to machine learning, AI etc. A cynical formulation would be to say that a research project is a machine for converting masses of data into complicated brightly coloured diagrams. The use of simple logical arguments to get real insights is becoming rarer. In my opintion this is not because such things are no longer possible or useful but decause they are no longer fashionable. Conversations with other older participants of the conference indicate to me that I am not alone in this opinion. Having got rid of some complaints let me now say something about some of the talks at the conference I liked best.

The first of these is a case where there were huge amounts of data involved and very complicated coloured pictures but there were also practical results which I found very impressive. The talk was by Bernd Bodenmiller from Zurich. In this work mass spectrometry techniques were used to produce very detailed pictures of the distribution of substances in slices of tumour tissue. I was surprised by one picture which showed the distribution of a conventional platinum-based chemotherapeutic agent within a tumour. While it was uniformly distributed through certain parts of the tumour it was more or less absent from others. Apparently cancer cells can develop methods to exclude this kind of drug from certain regions and thus survive. The one theme in the talk which caught my attention most was an application of this method to ovarial cancer. This is a very deadly cancer with frequent relapses after treatment. In the work reported on imaging techniques were use to distinguish different classes of patients and relate the differences between them to the rate of relapse. Beyond this predictions could be made which drugs might benefit which patients most. Patients were subjected to a kind of dual treatment strategy which would be illegal in Germany but which is fortunately legal in Switzerland. The idea is that on the one hand therapy decisions are considered in a conventional way and this information is given to the tumour board. Independently of this an analysis is done using the advanced imaging methods and this information also goes to the tumour board. These two sets of information are combined to make therapeutic decisions. In one case this method was applied to a patient already in palliative care, predicted to live for only a few more weeks. Five years later she is still alive and well. This is just one extreme case and the total sample size of patients is small. Nevertheless the preliminary conclusion is that this method leads to an large extension of the lifetime of the patients (I think a factor of four) in comparison to conventional approaches.

The second talk I want to mention is that of Becca Asquith. I had already heard a talk by her on a similar subject a couple of years ago and I wrote about it in a previous post. It has been observed that the KIRs an individual has can affect their ability to combat various infectious and autoimmune diseases, both positively and negatively, depending on the example. This is correlated to which MHC molecules the individual has. The subject of the talk was understanding the mechanisms behind these phenomena. One conclusion is that a determining factor is the typical lifetime of T cells. So how could KIRs modulate this lifetime? Two hypotheses are compared. One of these is that NK cells carrying the KIRs kill T cells, thus reducing their average lifetimes. Experiments were described which together with modelling, can decide between these two mechanisms. I find that this project was a beautiful combination of theory and experiment, exactly as I imagine such a project should ideally be. The whole thing, in particular the logical connections were very well described in the talk. At the end I asked the speaker why NK cells should kill T cells. Could this be of benefit to the organism or is it just a kind of collateral damage? My understanding of the answer, which I find plausible, is that any mechanism which can be used by the immune system to regulate its activity will be used.

The third talk was by Andreas Reichel, head of research at the company Bayer. He started off by mentioning a possible mechanism of action of a drug which is different to those commonly seen. This is to direct a certain protein to the proteasome so that it is destroyed. He then talked about the way in which candidate drugs are identified in the pre-clinical region. He mentioned a method in which a relatively simple ODE model can be used to obtain information. It can be used to find promising candidates. It can be used to suggest good doses for trials. (Sometimes increasing the dose produces no effect of the kind desired.) It can be used to choose optimal times for taking blood samples when testing candidates. Apparently it has been possible to convince decision makers that this theoretical work is something they can really profit from. For me this is a good example of how (relatively simple) mathematics can be used to make a significant contribution to a practical task such as drug discovery. If someone wants to develop models of this kind or apply them in an intelligent way then they need to things from analysis which I teach students on a day to day basis.

Catch bonds and T cell activation

November 17, 2025

A frequent approach to studying T cell activation is based on the idea that properties of the chemical binding of the T cell receptor to an antigen determine whether the cell is activated or not. This applies in particular to the work I have written about here and here. An alternative idea is that mechanical properties also play a role in this process. We have studied this in a new preprint with Yogesh Bali and Wolfgang Quapp. What happens when two molecules are chemically (non-covalently) bound and we apply a force which tends to pull them apart? A simple scenario would be to assume that the greater the force the greater the chance the bond will break. This is what is called a slip bond. However there is also another possibility. It may be that in a certain range the applied force actually stabilises the bond. This is what is called a catch bond. It is related to what is called a Chinese finger trap. This is a toy with the following property. If you push a finger into it it is difficult to get it out again and pulling harder only makes it worse. The way to escape is to push instead of pull.

In this work we study a model with two mechanical degrees of freedom which are the extension of the bond between T cell and the pMHC and the angle between the T cell receptor and the cell membrane. The dynamics is described by a potential depending on these two variables. Stable and unstable bonds correspond to minima and saddle points of this potential. The potential depends on a constant external force as a parameter and varying this parameter leads to bifurcations. In the paper this behaviour is studied while trying to incorporate as much experimental data as possible related to this system.

The mystery of thermodynamics

October 14, 2025

I have recently been reading the book ‘Foundations of Chemical Reaction Network Theory’ by Martin Feinberg. In the past I have spent a lot of time studying Feinberg’s classic lecture notes on Chemical Reaction Network Theory and I have read (parts of) many of his papers. Thus many things in the book are familiar to me. However there are also many things which are not and I feel that I am profiting a lot from reading it. One subject which plays a marginal role in most of Feinberg’s writings is thermodynamics and the related concept of detailed balancing. This is different in the book since there are two chapters (Chapter 13 and Chapter 14) devoted to these topics. On p. 274 we read ‘The mathematical foundations of thermodynamics remain somewhat murky, at least to me.’ My response to this statement is ‘me too’. In fact I have often experienced that mathematically inclined people say that they never understood thermodynamics. My own difficulties with the subject influenced my career. As a student I was irritated by the equation \frac{\partial V}{\partial T}\frac{\partial T}{\partial P} \frac{\partial P}{\partial V}=-1. When I asked the lecturer who was teaching us a course on thermodynamics he was not able to give me an explanation which I found satisfactory. As a schoolboy physics was the subject which interested me most. After my second year at university I had to decide between doing a degree in physics, a degree in mathematics or a joint degree in both. My decision for the second alternative was strongly influenced by that thermodynamics conundrum. Another experience which contributed to my decision was that at one time I happened to have two courses on the same topic, Fourier series, in the same term, one in physics and one in mathematics. The second was quite transparent to me, the first obscure. The fact that I had
a more positive experience with mathematics than with physics as a student probably had to do with the fact that the relative quality of the lecturers in mathematics was better. At the same time it has to with the nature of the subjects themselves. To come back to Feinberg’s book, on p. 281 he writes ‘When I was an undergraduate student, classical thermodynamics appeared to be a beautiful (and somewhat kabbalistic) subject, but its purpose was not clear. … I didn’t really understand what was happening.’ Feinberg and many other people seem to have made peace with thermodynamics although remaining with an uneasy feeling. This does not apply to me. Perhaps it is an aesthetic thing: Feinberg found the subject ‘beautiful’, even as a student, while I must say that I experienced it as ugly.

Thermodynamics is a part of physics which seems to be difficult to relate to rigorous mathematics. In this sense it bears a resemblance to the much more prominent example of quantum field theory. What does the word thermodynamics mean to me? I want to try to answer this question without reading what anyone else says about the subject. (I can do that later, if desired.) I start with an etymological approach. This indicates that the subject has to do with heat and the way that a system evolves in time. Another approach is a historical approach. I have the impression that a motivation for the subject was understanding the efficiency of steam engines. Yet another approach is to try to make contact to statistical mechanics. A gas is made up of an enormous number of molecules and it is impossible to keep track of them individually. Thus we pass to a statistical description. This involves some probability theory or possibly even quantum mechanics. Getting to thermodynamics involves discarding some information about the system and nevertheless ending up with a description which is to some extent self-contained.

Talk by Marisa Eisenberg at the SMB annual conference

July 16, 2025

At the moment I am attending the annual meeting of the Society for Mathematical Biology in Edmonton. Yesterday I heard a talk of Marisa Eisenberg which was partly about public health and partly about fundamental issues of mathematical modelling. I found both aspects of her talk interesting and also the connection between them. The speaker is Professor at the School of Public Health at the University of Michigan and Director of MICOM (Michigan Public Health Integrated Center for Outbreak Analytics and Modelling). She allowed herself an ironic comment on this acronym. In the context of this position she has had a lot of interactions with public health officials of the state of Michigan. At the beginning they did not appreciate what she had to tell them much. After a lot of talking this changed. Now they really take her input seriously. She emphasized that it can be useful to give these officials the kind of data they ask for even if it might be better to base a decision on other data. She stressed how important it is to make clear to people what modelling cannot do. Her activities have not been limited to the state level. She has also been working on a project to convince people of the value of measles vaccinations at a local level.

Turning to the more theoretical topics, one concept which was central in the talk was that of parameter identifiability. Very often in mathematical biology (and in epidemiology in particular) we have models consisting of a system of ordinary differential equations depending on parameters. Many of the parameters have not been measured directly. We can attempt to determine these parameters for a given application of the model by fitting to experimental data. However it can happen that in examples the parameters cannot even in principle be determined from measured data. Parameter identifiability can be defined as the absence of this problem. This is a subject which I have heard about in several talks in the past years, including some by Nicolette Meshkat. I found the topic interesting on an abstract mathematical level but never took it as seriously as I probably should have done from the point of view of its practical importance. This concept came up in many of the talks I have heard at this conference. Perhaps the time has come to end my relative negligence of this subject.

Sometimes even when the parameters are in principle determined by data the dependence of the data on the parameters is very weak so that identifying the parameters is impracticable. What can we do about this? One possibility mentioned in the talk is that the quantities of interest in the application are independent of the parameters which cannot be determined. Then we can simply ignore the problem without doing any damage. Another issue discussed in the talk was the comparison between mechanistic approaches and the use of artificial intelligence. A disadvantage of AI in epidemiology us that it needs a lot of training data and at the beginning of an epidemic not much data is available. An alternative put forward by the speaker is as follows. To understand a particular disease first produce a number of different mathematical models. Then use simulations of these models to produce data for AI. This could help to give AI a head start and partly overcome one of its biggest disadvantages in this context.

Backward bifurcations and multistationarity

May 20, 2025

In a previous post I wrote about backward bifurcations in the context of a model of hepatitis C. At that time I believed that this kind of bifurcation was not possible for that model. In recent work with Alexis Nangue we discovered that this belief was not correct. In a new preprint we prove that backward bifurcations do occur in that model. The existence of a backward bifurcation proves that there are parameters where the basic reproductive ratio {\cal R}_0 is less than one and there nevertheless exists a positive steady state. The bifurcation which occurs in the proof is a transcritical bifurcation. Intuitively this means that as a parameter is varied a steady state moves across another steady state which is fixed. There is a sign involved in a bifurcation of this type. One value of the sign corresponds to a backward bifurcation while the other corresponds to the forward bifurcations more familiar in epidemiology. The steady state which arises in a backward bifurcation is unstable. In simulations of models with backward bifurcations it is often seen that the steady state can be continued to larger values of the parameter until it reaches a point where a fold bifurcation takes place. There it changes from unstable to stable. This phenomenon cannot be captured by the analysis of the transcritical bifurcation since in general it occurs far away from the bifurcation point. We have now found an extension of that analysis which allows the fold bifurcation to be treated. The idea is to replace the transcritical bifurcation by one whose codimension is one higher. This is similar to passing from a fold to a cusp. The transcritical bifurcation could be given the intuitive name ‘moving hyperbolic steady state’ and by analogy with this we give the bifurcation studied in the paper the name ‘moving fold’. In this type of bifurcation two positive steady states arise, one stable and one unstable. This is a tool which may be used to prove the existence of more than one positive steady state in epidemiological models. The general theorem we prove applies to both the cases {\cal R}_0<1 and {\cal R}_0>1 but up to now we have only found interesting examples for the first of these cases. Our original model for hepatitis C is one such example. We were not able to decide whether in that example there exist parameters which would allow the second case of the theorem to be applied. We proved the relative statement that if parameters with this property exist there actually exist three positive steady states. An interesting side effect of this work is that we found a limiting case of that model (where two parameters are set to zero) where there can exist a continuum of steady states.

Mathematical models for T cell activation, part 2

April 12, 2025

In a previous post I discussed the topic of T cell activation and presented some results which Eduardo Sontag and I obtained on this subject. We studied a mathematical model of François et al. and found some interesting features of the solutions not previously known. We found that for some parameter values there exists more than one positive steady state. We also found that for some parameter values the response function describing the dependence of the activation on important parameters is non-monotone. It can happen that increasing the stimulation of the T cell decreases its level of activation. What we did not check was to what extent the parameter values for which these phenomena are observed are consistent with experimental measurements. Now in a paper of Yogesh Bali and myself we showed that they do occur for biologically reasonable parameter values. The paper also contains an extensive comparison of different models for the initial stages of T cell activation. In particular we give rigorous proofs of some properties of the solutions. We apply the Deficiency Zero Theorem in many cases. We also present simulations which allow comparisons of the predictions of different models with experimental results. A new discovery was that for one model the response function can have more than one maximum.

Matched asymptotic expansions, part 3

December 14, 2024

If the inner and outer solutions which we introduced previously are to define a good approximation then their difference should tend to zero as \epsilon tends to zero. Their derivatives should also do so. This idea is implemented by looking at some intermediate point which is not too close to the origin but not too far away so that it can be hoped that both approximations are valid at that point. In the book the choice x=\sqrt{\epsilon} is made and this is another case where intuition is used instead of a clear algorithm. This leads to the choice c=1 and d=0 and in that case the difference vanishes faster than any power of \epsilon as \epsilon\to 0. Since the equation is linear we can add the inner and outer solutions to get a new solution. Here the lowest order term is present in both summands and so we need to subtract it off. The result is y_u(x)=2x-1+e^{\frac{x}{\epsilon}}. It satisfies the boundary condition at x=0 exactly but it satisfies the boundary condition at x=1 only approximately, whereby the error is exponentially small.

What has been achieved here? The method of matched asymptotic expansions has been applied to identify a candidate approximation to the solution of a certain boundary value problem. Since in this example the solution can be computed explicitly it could be checked that this candidate actually does satisfy criteria which makes it a good solution. On the other hand if the exact solution is not used it gives us no independent argument that the candidate is a good approximation. An interesting question is how it is possible to obtain a rigorous proof that a function obtained by this method satisfies the conditions which it was designed to satisfy.

This problem can be related to the technique of geometric singular perturbation theory (GSPT). Let v=u'. Then the original ODE can be written in the form \epsilon v'=-v+2, u'=v. This is the canonical form of a system with two time scales in GSPT. We have a slow variable u and a fast variable v. The critical manifold is given by v=2. This problem is normally hyperbolic with an eigenvalue -1. Intuitively, the solution approaches the slow manifold very fast and then moves slowly along it. It is not possible to immediately apply the most basic techniques of GSPT to this problem because the solution on the slow manifold goes to infinity – it does not remain in a compact set.

Matched asymptotic expansions, part 2

December 13, 2024

In a previous post I wrote about matched asymptotic expansions, concentrating on a classical example, the Kaplun-Lagerstrom model. Now I want to go into some aspects of this subject more deeply. My main source is once again the book of Hastings and McLeod. In this post I discuss the ODE \epsilon y''+y'=2 on the interval [0,1] with boundary conditions y(0)=0 and y(1)=1. The equation contains a parameter \epsilon. We would like to understand the behaviour of solutions of this problem for \epsilon>0 small. The equation is linear and looks quite simple. That this problem is not straightforward can be seen if we set \epsilon=0. A naive hope would be that the solution for \epsilon small and positive can be approximated uniformly by a solution of the limiting equation with \epsilon=0. In fact the equation with \epsilon=0 has the general solution y=2x+c and this does not satisfy the desired boundary conditions for any value of c. This has to do with the fact that the equation for \epsilon>0, being second order, has a two-parameter family of solutions while that for \epsilon=0, being first order, has only a one-parameter family. This shows that this is what is called a singular perturbation problem, where the nature of the equation changes essentially as a certain parameter value is reached. To obtain more information we expand the solution in powers of \epsilon, y=\sum_{i=0}^\infty y_i\epsilon^i and substitute this into the equation. This means that we try the simplest expansion imaginable. Then y_0 satisfies the equation for \epsilon=0 and we are confronted with the problem which is already familiar. If we ignore this we find that it is consistent to set all y_i with i\ge 1 to zero. There is no clear preferred choice for the constant c. In fact the ODE and thus the boundary value problem can be solved explicitly. The strategy here is to ignore that fact as much as possible and gain information in other ways. The hope is that these other methods will apply in cases where no explicit solution is available. Looking at the explicit solution reveals that the choice c=-1 gives a function which solves the boundary condition y(1)=1 and is a good approximation to the true solution on most of the interval [0,1]. When x=0 is approached the exact solution swerves off from the approximate solution so as to head for the point corresponding to the boundary condition y(0)=0. The derivative becomes very large when \epsilon is very small. It looks like the solution could be got by piecing together an outer solution, the explicit one we already saw, with an inner solution defined in a small region near x=0.

How can we get information about the inner solution? One strategy is to use a magnifying glass to see more clearly what is happening in the small region of interest. This is done by introducing a new coordinate \zeta=\frac{x}{\epsilon}. This is a bit arbitrary. In principle we could have used a different power of \epsilon in the rescaling. It turns out that this choice is useful in this example. Call the transformed solution z(\zeta). It satisfies \frac{d^2 z}{d\zeta^2}+\frac{dz}{d\zeta}=2\epsilon. We can now do an expansion of this equation in powers of \epsilon. Substituting this in the original ODE gives a sequence of equations for the expansion coefficients z_i. The case i=0 can be solved to give z_0(\zeta)=c(e^{-\zeta}-1). This is the general solution satisfying z_0(0)=0. In the case i=1 the general solution with z_1(0)=0 is of the form z_1(\zeta)=2\zeta+d(e^{-\zeta}-1). It is consistent to set z_i=0 for i\ge 2 and we do so. Thus we have two free constants c and d. We now want to combine these inner and outer solutions to approximate the solution of the original problem as well as possible. We try to match them as well as possible by choosing the parameters c and d. An account of how this is done will be postponed to a later post.

Bessel functions

November 24, 2024

When I am lecturing on elementary topics related to ordinary differential equations I always emphasize that most ODE cannot be solved explicitly. What is an ‘explicit’ or ‘exact’ solution? It means that the solution can be written as an algebraic expression containing well-known (‘elementary’) functions such as trigonometric functions, exponential, logarithm etc. So what can we do when no explicit solution exists? Option 1: give up. Option 2: try to calculate the solution numerically on a computer. Option 3: apply heuristic techniques such as assuming that some quantities are ‘small’ and then setting them to zero to get a simpler equation. Note that the result of Options 2. and 3. are functions which are supposed to be approximations to the solution in some (often unspecified) sense. Option 4 (my favourite) is to study the qualitative properties of solutions. This does include the possibility of assuming that some parameters are small in a precise sense to get a simplified equation and then proving that under certain circumstances solutions of the limiting equation approximate solutions of the original equation in a precise sense. To avoid any misunderstanding I regard all these options (except Option 1) as often being valuable and I believe that combining them gives the best result. We can think of an ‘elementary function’ as being the solution of some particularly simple ODE. Since these solutions have often been studied it makes sense to give them a name. They can then be studied in cases where they arise by exploiting the knowledge obtained by previous generations of mathematicians.

There is always the option of extending the class of elementary functions by defining a class of functions which are the solution of a certain ODE. Here I want to look at an example of this, the Bessel functions. They are the solutions of the homogeneous linear equation x^2\frac{d^2y}{dx^2}+x\frac{dy}{dx}+(x^2-\alpha^2)y=0. Here \alpha is an arbitrary complex number and we also think of x and y as being complex. These functions are probably more familiar to physicists than to mathematicians since they arise when looking for solutions of certain PDE important in physics (e.g. the Laplace equation) which inherit a symmetry property of the equation (e.g. spherical or cylindrical symmetry). Evidently the space of solutions of the Bessel equation is a two-dimensional vector space. The equation is symmetric under a change of sign of x and so we may as well restrict consideration to the region x\ge 0. It is also symmetric under a change of sign of \alpha. The solutions behave in an oscillatory manner. For large values of x they look like x^\alpha times an oscillatory function of order one. In particular, if \alpha<0 they have a power-law decay. For small values of x some are unbounded, of order x^\alpha, while others are bounded. There are a lot of notations related to Bessel functions and I refer to the corresponding Wikipedia article for an overview. I will mainly use the notation J_\alpha (x). It should be noted that this notation includes a choice of the sign of \alpha, so that J_\alpha and J_{-\alpha} are often different, although they solve the same equation. It also includes a choice of the normalization of the function. The asymptotics for x\to 0 depend on \alpha. For \alpha>0 or \alpha an integer J_\alpha is bounded near the origin. For negative non-integer \alpha it is unbounded. For non-integer \alpha the functions J_\alpha and J_{-\alpha} are linearly independent. This means that if \alpha>0 in this case J_\alpha is up to normalization the only bounded solution. For \alpha an integer we have J_{-\alpha}=(-1)^\alpha J_\alpha and so the J_\alpha, which are sometimes called Bessel functions of the first kind, only provide a one-dimensional solution space in that case. To get the full solution space it is necessary to introduce another linearly independent solution, the Bessel function of the second kind Y_\alpha. The function Y_\alpha can be defined algebraically in terms of J_{\alpha} for \alpha not an integer and then by taking a limit of this for \alpha an integer. The Y_\alpha are singular at x=0. They satisfy Y_\alpha=(-1)^\alpha Y_{-\alpha} when \alpha is an integer.

There is a series expansion for the J_\alpha which may be found in the Wikipedia article. Here I just note that up to a multiplicative constant the leading order term near x=0 is x^\alpha. The higher order terms are proportional to x^{\alpha+2m} for integer m. Thus the boundedness properties of these functions near the origin correspond to those of the complex powers. For \alpha a positive integer the leading order term in Y_\alpha is a multiple of x^{-\alpha} and there is a higher order term proportional to x^{\alpha}\log x. For \alpha=0 the leading term is a logarithm.

What does the deficiency depend on?

August 12, 2024

I have spent a lot of time and energy thinking about chemical reaction networks and, in particular, about the notion of the deficiency of such networks. This is a non-negative number which can be assigned to a network and it gives a certain measure of the complexity of the network. Now I realised something which I should have realised long ago. I am sure that it is widely known and plays a role, at least implicitly, in many papers I have read and talks I have heard. As one recent example of a relevant text I mention a paper of Craciun et al., Math. Biosci. 342, 108720. Given a reaction network describing a system of chemical reactions there is a standard way of obtaining a corresponding system of ODE, assuming mass action kinetics. Given the network we can also compute the deficiency. The point I had not fully understood is that the deficiency cannot be computed from the ODE system. It is a property of the network and more than one network can give rise to the same system of ODE. In particular, two networks with different deficiencies can give rise to the same system of ODE. Thus it makes no sense to talk about the deficiency of the system of differential equations.

A situation in which this can play a role is the following. The ODE corresponding to a reaction network can have linear conserved quantities. In the usual language of chemical reaction network theory (CRNT) the stoichiometric subspace is a proper subspace. It is then often convenient to eliminate some of the variables so as to get an effective system of lower dimension, thus obtaining a simplification of the system. It should however be noted that what is being reduced here is the system of ODE, not the network. Thus I cannot speak about the deficiency of the reduced network. I must first find a network which gives rise to the reduced system and then calculate its deficiency. In general there will be a choice involved in defining the network and the deficiency which comes out may depend on that choice.


Design a site like this with WordPress.com
Get started