# Intersection Between Stable and Saddle-Point Solution -Stability Analysis

I have a free energy function:

$$G(N_b, l_b)= -N_b E_b + \frac{1}{2} N_b \kappa_b (l_b - 1)^2 - F(l_b - 1) + \frac{A}{2} k_g (u_g - l_b)^2 + (N_t - N_b) \text{Log}\big(\frac{N_t - N_b}{A}\big) + N_b \text{Log}\big(\frac{N_b}{A}\big)$$

Physically, this represents the fee energy function of two substrates conjugated by $$N_b$$ reversible bonds that are being pulled away by a constant force $$F$$. Minimizing the free energy $$G$$ with respect to the state variables $$N_b$$ and $$l_b$$ ($$\partial G/\partial N_b=0$$ and $$\partial G/\partial l_b=0$$) yields the equilibrium equations for the system. I have solved these equations numerically, using FindRoot. For a given force $$F$$, I obtain two roots, one is the stable solution (solid-line branch) and the other one is the saddle-point solution (dashed-line branch) as shown in the figure below. I used the second derivative test (checking the signs of Hessian eigenvalues) to make sure these roots are indeed the local minima and saddle-point in the free energy landscape.

You can see that there is a critical force $$F_{cr}$$ where these two branches intersect; this is the force at which the assembly snaps and the two substrates completely separate. What I want to obtain is an expression for $$F_{cr}$$ as a function of one of the input variables, such as $$\kappa_b$$. More specifically, my question is that mathematically what happens at the intersection. Of course, finding an analytical expression of $$F_{cr}$$ as a function of $$\kappa_b$$ would be ideal but using a numerical method to get $$F_{cr}$$ vs $$\kappa_b$$ curve would be great too. The curves shown above are constructed using $$N_t=50$$, $$E_b=10$$, $$k_g=2.43902\times10^{-6}$$, $$u_g=3$$, and $$A=100$$.

The Mathematica code to obtain the graphs shown above is simple. I used FindRoot to solve the equilibrium equations for a given F and used the results as the initial values for the next FindRoot at a slightly larger F.


Clear["Global*"]

G[Nb_, lb_] := -Nb Eb + 1/2 Nb kb (lb - 1)^2 - F (lb - 1) +
A/2 kg (ug - lb)^2 + (Nt - Nb) Log[(Nt - Nb)/A] + (Nb) Log[Nb/A]

D[G[Nb, lb], Nb]

-Eb + 1/2 kb (-1 + lb)^2 + Log[Nb/A] - Log[(-Nb + Nt)/A]

D[G[Nb, lb], lb]

-F + kb (-1 + lb) Nb - A kg (-lb + ug)

A = 10^2; Nt = 50; kb = 24.39; kg = 2.43*^-6; Eb = 10; ug = 3;

F = 1;

Ff[1] = F;

Clear[x,z]

FindRoot[{-Eb + 1/2 kb (-1 + lb)^2 + Log[Nb/A] - Log[(-Nb + Nt)/A] ==
0, -F + kb (-1 + lb) Nb - A kg (-lb + ug) == 0}, {{lb,
1.0008}, {Nb, 49.99}}]

{lb -> 1.00082, Nb -> 49.9977}

x = lb /. %; z = Nb /. %;

lb[1] = x; Nb[1] = z;

F = 10;

Ff[2] = F;

Clear[x,z]

FindRoot[{-Eb + 1/2 kb (-1 + lb)^2 + Log[Nb/A] - Log[(-Nb + Nt)/A] ==
0, -F + kb (-1 + lb) Nb - A kg (-lb + ug) == 0}, {{lb,
1.0008204438072899}, {Nb, 49.99773008793257}}]

{lb -> 1.0082, Nb -> 49.9977}

x = lb /. %; z = Nb /. %;

lb[2] = x; Nb[2] = z;


The same equations are used to obtain the saddle-point solutions. The only difference is the initial values used for FindRoot. I understand that the manual updating of the initial values for FindRoot is rather inefficient but I couldn't figure a better way to do it. At the end, I plotted the collected lb[i] and Nb[i] values to get the graphs shown above.

• Please include some Mathematica code, at the very minimum for your expressions. Commented Aug 20, 2023 at 13:32
• I added a piece of my code. Hopefully that is enough. Commented Aug 20, 2023 at 14:11
• What is A1 in your code definition of G? I would also suggest commenting out code output (*like this*) in your code block so that it doesn't get evaluated as input when copy-pasting. It also makes it easier to read.
– ydd
Commented Aug 20, 2023 at 16:14
• I don't understand how you get the dashed lines. Could you include the code that generates one of the curves? Commented Aug 20, 2023 at 18:50

Your original statement I seems to have solved analytically, but the solution is so large that does not admit any analysis. I propose that you limit the model. You did not explain the physical sense of the terms 1/2 Nb kb (lb - 1)^2 and A/2 kg (ug - lb)^2 but it seems that they play comparable roles: the energies of the bonds in the spring approximation.

The first thing I suggest: let us for the first study omit one of them, say, A/2 kg (ug - lb)^2 The second suggestion is to reduce the number of parameters by the rescaling:

rule = {Nb -> n*A, Nt -> nt*A, F -> f*A, Eb -> e*kb, lb -> l - 1};


Then introducing your free energy without the discussed term:

G[Nb_, lb_] := -Nb Eb + 1/2 Nb kb (lb - 1)^2 - F (lb - 1) +
(Nt - Nb) Log[(Nt - Nb)/A] + (Nb) Log[Nb/A]


and rescaling one gets the following:

g = G[Nb, lb]/ A   /. rule // Simplify

(* 2 f - f l + 2 kb n - e kb n - 2 kb l n + 1/2 kb l^2 n +
n Log[n] + (-n + nt) Log[-n + nt]  *)


Here are equations of state:

eq1 = D[g, n] == 0
eq2 = D[g, l] == 0

(*  2 kb - e kb - 2 kb l + (kb l^2)/2 + Log[n] - Log[-n + nt] == 0

-f - 2 kb n + kb l n == 0  *)


The solution for l for future use:

slL = Solve[eq2, l][[1, 1]]

(*  l -> (f + 2 kb n)/(kb n)  *)


Let us substitute it to the first equation:

eq1A = eq1 /. slL // Simplify

(* f^2/(2 kb n^2) + Log[n] == e kb + Log[-n + nt] *)


and solve it with respect to the rescaled force, f:

Solve[eq1A, f][[2, 1]]

(* f -> Sqrt[2] Sqrt[kb] n Sqrt[e kb - Log[n] + Log[-n + nt]] *)


Let us have a look at the obtained solution at arbitrary values of all parameters. Namely I put them all equal to 1:

Plot[Sqrt[2] Sqrt[kb] n Sqrt[
e kb - Log[n] + Log[-n + nt]] /. {kb -> 1, nt -> 1, e -> 1}, {n, 0,
1}, AxesLabel -> {Style["n", 16, Black, Italic],
Style["f", 16, Black, Italic]}]


The maximum of this curve is what you are looking for. It is given by the following equation:

eqMax = Numerator[
D[ n Sqrt[e kb - Log[n] + Log[-n + nt]], n] // Simplify] == 0

(*  -1 + 2 e kb + n/(n - nt) - 2 Log[n] + 2 Log[-n + nt] == 0  *)


Expressing n/(nt - n)==x as follows

slN = Solve[n/(nt - n) == x, n][[1, 1]]

(*  n -> (nt x)/(1 + x)  *)


and substituting this into an equation one gets:

eqMax1 = eqMax /. slN // Simplify

(*  2 e kb == 1 + x + Log[x^2]  *)


This equation one can solve analytically:

Solve[eqMax1, x] // Quiet

(*  {{x -> 2 ProductLog[-(1/2) Sqrt[E^(-1 + 2 e kb)]]}, {x ->
2 ProductLog[1/2 Sqrt[E^(-1 + 2 e kb)]]}}   *)


whence

Solve[n/(nt - n) == 2 ProductLog[1/2 Sqrt[E^(-1 + 2 e kb)]], n]

(*  {{n -> (2 nt ProductLog[1/2 Sqrt[E^(-1 + 2 e kb)]])/(
1 + 2 ProductLog[1/2 Sqrt[E^(-1 + 2 e kb)]])}}   *)


Is the exact analytical solution of the equation.

Further, substituting this solution into the equations eq1 and eq2 one finds:

((2 kb - e kb - 2 kb l + (kb l^2)/2 + Log[n] - Log[-n + nt] /. slN1 //
FullSimplify) /. -Log[a_] + Log[b_] -> Log[b/a] //
FullSimplify) == 0

(*  1/2 kb (-2 e + (-2 + l)^2) +
Log[2 ProductLog[1/2 Sqrt[E^(-1 + 2 e kb)]]] == 0  *)


and

-f - 2 kb n + kb l n == 0 /. slN1 // FullSimplify

(*  f + kb (-2 + l) nt (-1 + 1/(
1 + 2 ProductLog[1/2 Sqrt[E^(-1 + 2 e kb)]])) == 0  *)
`

yields you two conditions at which the critical point is achieved. One can still play with these equations and try to build the manifold on which this happens.
As much as I understand, mathematically, on this manifold, one of two minimums of the free energy fuses with the maximum situated between them. As a result, only one free energy minimum survives.

• Thank you, Alexei. It's quite expected, considering your profound understanding of cell mechanics, assuming you're indeed the Alexei Boulbitch I'm referring to. Commented Aug 20, 2023 at 20:08
• You are welcome. Commented Aug 21, 2023 at 9:09