UNIVERSITY OF TWENTE.

MATHEMATICAL METHODS ASSIGNEMENT 1

BY - ERMIAS ATNAFU
Student number – 3191656

Date of Submission - December 3, 2024

1
 Wilson-Cowan: A model of activity of neural populations

1. Write a matlab function WC.m so that the call f=WC(t,X) is the evaluation of (1) and f = (E0 , I0 ) is a column
vector.

Matlab Code:

function f = WC(t, X)

% Parameters
c1 = 20;
c2 = 11.5;
c3 = 15;
c4 = 1;
P = 2;
Q = 1;
bE = 1.3;
bI = 2;
thetaE = 4;
thetaI = 3.7;

% Extract E and I from X

E = X(1);
I = X(2);

% Define the sigmoid functions

SE = @(u) 1 / (1 + exp(-bE * (u - thetaE))) - 1 / (1 + exp(bE * thetaE));
SI = @(u) 1 / (1 + exp(-bI * (u - thetaI))) - 1 / (1 + exp(bI * thetaI));

% Evaluate the equations

E0 = -E + (1 - E) * SE(c1 * E - c2 * I + P);
I0 = -I + (1 - I) * SI(c3 * E - c4 * I + Q);

% Return the result as a column vector

f = [E0; I0];
end

2. Write a matlab function JACWC.m so that the call A=JACWC(t,X) returns the Jacobian matrix fX. We expect you
to use the result from Tutorial 1.

As per Tutorial 1 the MATLAB CODE will be :

%% Exercise 2

function A = JACWC(t, X, h)

% JACWC: Compute the numerical Jacobian matrix for the Wilson-Cowan system

% Inputs:

% t - Time variable (not used but required for compatibility)

2
 % X - State vector [E; I] where E and I are the excitatory and inhibitory populations

% h - Small perturbation for numerical differentiation

% Output:

% A - Numerical Jacobian matrix

n = length(X); % Dimension of the state vector

A = zeros(n, n); % Initialize the Jacobian matrix

for i = 1:n

v = zeros(n, 1); % Perturbation vector

v(i) = 1; % Perturb one variable at a time

% Evaluate the perturbed values of the function

X_plus = WC(t, X + v * h);

X_minus = WC(t, X - v * h);

% Compute the numerical derivative (finite differences)

A(:, i) = (X_plus - X_minus) / (2 * h);

end

end

%% Example Usage

% Define the Wilson-Cowan model function (already created above as WC.m)

% Ensure WC.m is in the same directory or MATLAB path.

% Set an example equilibrium point (you may calculate or assume this)

eq_points = [0.2, 0.1]; % Example equilibrium point [E*; I*]

X = [eq_points(1), eq_points(2)]; % State vector at equilibrium

% Set a small perturbation value

h = 1e-3; % Small step size for numerical differentiation

% Evaluate the Jacobian matrix at the equilibrium point

t = 0; % Current time (irrelevant here but necessary for WC's input)

A = JACWC(t, X, h);

3
 % Display the numerical Jacobian matrix

disp('Numerical Jacobian:');

disp(A);

%% Analytical Jacobian (for verification)

% Define partial derivatives analytically:

% J = [dE/dE, dE/dI;

% dI/dE, dI/dI]

theta_e = 4; be = 1.3;

theta_i = 3.7; bi = 2;

P = 2; Q = 1;

c1 = 20; c2 = 11.5;

c3 = 15; c4 = 1;

% Define helper functions

Se_prime = @(u) be * exp(-be * (u - theta_e)) / (1 + exp(-be * (u - theta_e)))^2;

Si_prime = @(u) bi * exp(-bi * (u - theta_i)) / (1 + exp(-bi * (u - theta_i)))^2;

% Compute partial derivatives

u_e = c1 * eq_points(1) - c2 * eq_points(2) + P;

u_i = c3 * eq_points(1) - c4 * eq_points(2) + Q;

dEdE = -1 + (1 - eq_points(1)) * Se_prime(u_e) * c1 - Se(u_e);

dEdI = -(1 - eq_points(1)) * Se_prime(u_e) * c2;

dIdE = (1 - eq_points(2)) * Si_prime(u_i) * c3;

dIdI = -1 - (1 - eq_points(2)) * Si_prime(u_i) * c4;

% Construct the analytical Jacobian

Jac = [dEdE, dEdI;

dIdE, dIdI];

% Display the analytical Jacobian

disp('Analytical Jacobian:');
disp(Jac);

4
 3. The equations dE/dt = 0 and dI/dt = 0 define the nullclines. On the nullclines, solution curves turn direction, i.e.
up/right ↔ down/left. Since E0 = 0 involves only one term with I, we can write I as a function of E. Plot the
resulting curves. (1.5pt) Alternative to finding an analytical expression for the nullclines, you may consider a
root finding method. In this case, for every fixed value of E, find I such that E0 = 0, and do this for a suitable
range of values for E.

MATLAB CODE :

% Parameters
b_e = 1.3; b_i = 2; % Sigmoid slopes
theta_e = 4; theta_i = 3.7; % Sigmoid thresholds
c1 = 20; c2 = 11.5; % Coupling constants for E
c3 = 15; c4 = 1; % Coupling constants for I
P = 2; Q = 1; % Inputs

% Pre-computed constants for sigmoid

delta_e = theta_e;
delta_i = theta_i;
exp_e = 1 / (1 + exp(b_e * delta_e));
exp_i = 1 / (1 + exp(b_i * delta_i));

% Define ranges for E and I

E = 0.01:0.01:0.49; % Excitatory activity range
I = 0.01:0.01:0.5; % Inhibitory activity range

% Define nullclines
I_E = @(E) (log(1 ./ (E ./ (1 - E) + exp_e) - 1) ./ b_e + c1 * E + P - delta_e) / c2;
E_I = @(I) (-log(1 ./ (I ./ (1 - I) + exp_i) - 1) ./ b_i + c4 * I - Q + delta_i) / c3;

% Compute nullclines
I_e = I_E(E); % Excitatory nullcline
E_i = E_I(I); % Inhibitory nullcline

% Plotting the nullclines

figure(1);
plot(E, I_e, 'b-', 'LineWidth', 2); % Excitatory nullcline
hold on;
plot(E_i, I, 'r--', 'LineWidth', 2); % Inhibitory nullcline
grid on;

% Labels and title

xlabel('Excitatory Activity (E)', 'FontSize', 12);
ylabel('Inhibitory Activity (I)', 'FontSize', 12);
title('Nullclines of the Wilson-Cowan Model', 'FontSize', 14);

% Legend
legend('Excitatory Nullcline', 'Inhibitory Nullcline', 'Location', 'southeast');
hold off;

5
 4. Find all equilibria of the Wilson-Cowan model for these parameter values numerically.(1pt)

MATLAB CODE :

% Parameters
b_e = 1.3; b_i = 2; % Sigmoid slopes
theta_e = 4; theta_i = 3.7; % Sigmoid thresholds
c1 = 20; c2 = 11.5; % Coupling constants for E
c3 = 15; c4 = 1; % Coupling constants for I
P = 2; Q = 1; % Inputs

% Pre-computed constants for sigmoid

delta_e = theta_e;
delta_i = theta_i;
exp_e = 1 / (1 + exp(b_e * delta_e));
exp_i = 1 / (1 + exp(b_i * delta_i));

% Define the system of equations

nullcline_system = @(X) [
(log(1 / (X(1) / (1 - X(1)) + exp_e) - 1) / b_e + c1 * X(1) + P - delta_e) / c2 - X(2);
(-log(1 / (X(2) / (1 - X(2)) + exp_i) - 1) / b_i + c4 * X(2) - Q + delta_i) / c3 - X(1)
];

% Initial guesses for the equilibrium points

initial_guesses = [
0.1, 0.1;
0.4, 0.4;
0.2, 0.3;
0.3, 0.2
6
];

% Find the equilibrium points using fsolve

options = optimoptions('fsolve', 'Display', 'off');
equilibrium_points = [];

for i = 1:size(initial_guesses, 1)
[eq_point, fval, exitflag] = fsolve(nullcline_system, initial_guesses(i, :), options);
if exitflag > 0
equilibrium_points = [equilibrium_points; eq_point];
end
end

% Remove duplicate points

% Define the sigmoid functions

S_e = @(E, I) 1 / (1 + exp(-b_e * (c1 * E - c2 * I + P - delta_e))) - 1 / (1 + exp(b_e *
delta_e));
S_i = @(E, I) 1 / (1 + exp(-b_i * (c3 * E - c4 * I + Q - delta_i))) - 1 / (1 + exp(b_i *
delta_i));

% Define the derivatives

E_prima = @(E, I) -E + (1 - E) * S_e(E, I);
I_prima = @(E, I) -I + (1 - I) * S_i(E, I);

dSedE = @(E, I) b_e * S_e(E, I) * (1 - S_e(E, I)) * c1;

dSedI = @(E, I) b_e * S_e(E, I) * (1 - S_e(E, I)) * -c2;
dSidE = @(E, I) b_i * S_i(E, I) * (1 - S_i(E, I)) * c3;
dSidI = @(E, I) b_i * S_i(E, I) * (1 - S_i(E, I)) * -c4;

dEdE = @(E, I) -1 - S_e(E, I) + (1 - E) * dSedE(E, I);

dEdI = @(E, I) (1 - E) * dSedI(E, I);
dIdE = @(E, I) (1 - I) * dSidE(E, I);
dIdI = @(E, I) -1 - S_i(I, I) + (1 - I) * dSidI(E, I);

% Initialize equilibrium points

eq_points = [];

% Define possible equilibrium points (initial guesses)

possible_eq_points = [
0.1, 0.1;
0.4, 0.4;
0.2, 0.3;
0.3, 0.2
];

for i = 1:size(possible_eq_points, 1)
X = possible_eq_points(i, :);
count = 0;

while true
% Compute function and Jacobian
F = [E_prima(X(1), X(2)); I_prima(X(1), X(2))];
Jac = [dEdE(X(1), X(2)), dEdI(X(1), X(2));
dIdE(X(1), X(2)), dIdI(X(1), X(2))];

% Solve linear system for correction

h = -Jac \ F;
X = X + h;

% Check convergence
7
 if norm(F) < tol
is_new_point = true;
if ~isempty(eq_points)
for k = 1:size(eq_points, 1)
if norm([X(1) X(2)] - eq_points(k, :)) < tol
is_new_point = false;
break;
end
end
end
if is_new_point
eq_points = [eq_points; [X(1) X(2)]];
end
break;
end
if count > 15000
break;
end
count = count + 1;
end
end

disp(eq_points)

Result: [0.1541 0.1921; 0.4727 0.4997 ; 0.4346 0.4995] these are the equilibrium points.

5. • Determine the linearization for each fixed point. Classify the fixed points based on the eigenvalues. Make local
phase portraits of the linearization x 0 = Df(E0, I0)x in separate figures.

MATLAB CODE:
eig_vals = [];

for i = 1:length(eq_points)

% Extract fixed point coordinates

X = [eq_points(i, 1); eq_points(i, 2)];

% Compute Jacobian matrix

Jac = [dEdE(X(1), X(2)), dEdI(X(1), X(2));

dIdE(X(1), X(2)), dIdI(X(1), X(2))];

a = Jac(1, 1);

b = Jac(1, 2);

c = Jac(2, 1);

d = Jac(2, 2);

% Compute eigenvalues

eigenvalue_1 = ((a + d) + sqrt((a + d)^2 - 4 * (a * d - b * c))) / 2;

eigenvalue_2 = ((a + d) - sqrt((a + d)^2 - 4 * (a * d - b * c))) / 2;

8
 eig_vals = [eig_vals; [eigenvalue_1, eigenvalue_2]];

% Classify fixed point

if imag(eigenvalue_1) == 0 && imag(eigenvalue_2) == 0

if real(eigenvalue_1) < 0 && real(eigenvalue_2) < 0

fprintf('The fixed point [%f, %f] is a Stable Node.\n', X(1), X(2));

elseif real(eigenvalue_1) > 0 && real(eigenvalue_2) > 0

fprintf('The fixed point [%f, %f] is an Unstable Node.\n', X(1), X(2));

else

fprintf('The fixed point [%f, %f] is a Saddle Point.\n', X(1), X(2));

end

elseif imag(eigenvalue_1) ~= 0

if real(eigenvalue_1) < 0 && real(eigenvalue_2) < 0

fprintf('The fixed point [%f, %f] is a Stable Spiral Focus.\n', X(1), X(2));

elseif real(eigenvalue_1) > 0 && real(eigenvalue_2) > 0

fprintf('The fixed point [%f, %f] is an Unstable Spiral Focus.\n', X(1), X(2));

else

fprintf('The fixed point [%f, %f] is a Center.\n', X(1), X(2));

end

end

end

Result:

• Fixed Point [0.1541, 0.1921]: Unstable Spiral Focus.

• Fixed Point [0.4727, 0.4997]: Stable Node.

• Fixed Point [0.4346, 0.4995]: Saddle Point.

6. Make a global phase portrait for a relevant part of the unit square 0 ≤ E, I ≤ 1. Your figure should contain the
nullclines, the equilibria and a few solutions (orbits). Create these solutions by simulating the model with the
Euler-Forward method and use at most 10 initial conditions. For the equilibria with real eigenvalues, also plot
the lines given by X(s) = x0+sv with v an eigenvector. Provide a short description of this figure and comment on
the attractors, i.e., where do the orbits go to?

Matlab code:

9
% Parameters
dt = 0.01; % Time step for Euler - Forward method
t_end = 10; % Simulation duration
t_steps = t_end / dt;

% Simulate Orbits
initial_conditions = [0.2, 0.2; 0.8, 0.2; 0.2, 0.8; 0.8, 0.8;
0.5, 0.5; 0.1, 0.9; 0.9, 0.1; 0.4, 0.6;
0.6, 0.4; 0.5, 0.3; 0.3, 0.5; 0.2, 0.1;
0.1, 0.2; 0.9, 0.5; 0.5, 0.9; 0.7, 0.8;
0.8, 0.7];
orbits = cell(size(initial_conditions, 1), 1);
for i = 1:size(initial_conditions, 1)
E = zeros(1, t_steps);
I = zeros(1, t_steps);
E(1) = initial_conditions(i, 1);
I(1) = initial_conditions(i, 2);
for t = 1:t_steps-1
dE = E_prima(E(t), I(t));
dI = I_prima(E(t), I(t));
E(t+1) = E(t) + dE * dt;
I(t+1) = I(t) + dI * dt;
end
orbits{i} = [E; I];
end

% Define ranges for E and I

E_range = 0.01:0.01:0.49;
I_range = 0.01:0.01:0.5;

% Compute nullclines
I_E = @(E) (log(1 ./ (E ./ (1 - E) + exp_e) - 1) ./ b_e + c1 * E + P - delta_e) / c2;
E_I = @(I) (-log(1 ./ (I ./ (1 - I) + exp_i) - 1) ./ b_i + c4 * I - Q + delta_i) / c3;
E_zero = [E_range', I_E(E_range)']; % Excitatory nullcline
I_zero = [E_I(I_range)', I_range']; % Inhibitory nullcline

% Create grid for quiver plot

[E_grid, I_grid] = meshgrid(0:0.05:1, 0:0.05:1);
E_prime = zeros(size(E_grid));
I_prime = zeros(size(I_grid));
for i = 1:numel(E_grid)
dE = E_prima(E_grid(i), I_grid(i));
dI = I_prima(E_grid(i), I_grid(i));
norm_factor = sqrt(dE^2 + dI^2);
if norm_factor ~= 0
distance_to_eq = min(sqrt((E_grid(i) - eq_points(:, 1)).^2 + (I_grid(i) - eq_points(:,
2)).^2));
scaling_factor = 1 / (1 + distance_to_eq);
E_prime(i) = (dE / norm_factor) * scaling_factor * 0.2; % Shorten vectors
I_prime(i) = (dI / norm_factor) * scaling_factor * 0.2; % Shorten vectors
end
end

% Plotting
figure;
hold on;
plot(E_zero(:, 1), E_zero(:, 2), 'b-', 'LineWidth', 2, 'DisplayName', 'Excitatory Nullcline');
plot(I_zero(:, 1), I_zero(:, 2), 'r--', 'LineWidth', 2, 'DisplayName', 'Inhibitory Nullcline');

% Plot Equilibria

10
scatter(eq_points(:, 1), eq_points(:, 2), 50, 'k', 'DisplayName', 'Equilibria');

% Plot Orbits
for i = 1:length(orbits)
orbit = orbits{i};
plot(orbit(1, :), orbit(2, :), 'k', 'HandleVisibility', 'off');
end

% Plot Quiver
quiver(E_grid, I_grid, E_prime, I_prime, 'Color', [210/255, 180/255, 140/255], 'LineWidth', 1,
'MaxHeadSize', 0.5, 'DisplayName', 'Vector Field');

% Final touches
xlabel('E');
ylabel('I');
xlim([0 1]);
ylim([0 1]);
title('Global Phase Portrait');
legend show;
grid on;
hold off;

11
 Figure Description

The output plot, titled final plot.jpg, presents a global phase portrait for the Wilson-Cowan model within the
bounds of (0 \leq E, I \leq 1). The plot includes the following elements:

 Nullclines: The nullclines for the excitatory and inhibitory populations are depicted, indicating where the
derivatives (E') and (I') are zero. These curves are shown between 0 and 0.5, as only the non-imaginary values
are plotted.

 Equilibria: The fixed points identified in Exercise 4 are marked as black circles, representing the system's
equilibrium points.

 Eigenvector Lines: For fixed points with real eigenvalues, lines representing the eigenvectors are drawn,
illustrating the directions of the local phase space near each equilibrium.

 Orbits: The trajectories for various initial conditions are plotted, showing the system's time evolution from
different starting points.

3. Write a report about the methods you have used and the results obtained. Motivate your choice for the algorithms,
step sizes, tolerances and other relevant parameters. Comment in detail on the results obtained. Do they make sense,
and are they accurate? Also, hand in the working code (in electronic form) that you have used to obtain your results.
We should not need to debug your code. There should be a single main-script that reproduces all your results. Call it
main with your initials as a suffix, e.g., I would write main HGEM. The code should be readable, i.e. contain
documentation. It should be clear what you do at every step by choosing self-explanatory names for your variables and
referencing algorithms in the book or lectures. If parts of your solution use the symbolic toolbox or commands such as
fsolve, eig or equivalent performing the numerical work in a ”black box”, this part will not be given credit. While we
check the code, ideally, the report should be self-contained.

Report on Phase Portraits of the Model

1. Introduction
This report presents the methods and results for simulating and visualizing the behavior of a two-variable dynamical
system using the Euler method for numerical integration. The system is governed by two functions, representing the
time derivatives of excitatory and inhibitory variables. We aim to explore the system's phase portrait by computing the
nullclines, equilibrium points, local phase portraits, and global phase portraits.
2. Methodology
2.1. WC Functions
The system of ordinary differential equations (ODEs) is defined by two functions, E_prima(E, I) and I_prima(E, I),
representing the rate of change of the excitatory (E) and inhibitory (I) variables. These functions are derived based on
the biological model that describes the interaction between excitatory and inhibitory neurons. We use these functions
to simulate the trajectories of the system over time using numerical methods.
2.2. Numerical Integration
To simulate the orbits of the system, we used the Euler method for numerical integration. The Euler method is simple
and efficient for solving ODEs but may accumulate errors over time, especially for stiff systems. However, given the
relatively small time step (dt = 0.01), the Euler method provides reasonable accuracy for the given problem. The time
step dt was chosen small enough to ensure smooth trajectories while balancing computational efficiency.
The duration of the simulation is t_end = 10, and the number of time steps is determined by dividing t_end by dt. The
initial conditions for the excitatory and inhibitory variables are chosen arbitrarily to cover a variety of initial states.
2.3. Nullclines
Nullclines are curves in the phase plane where the derivative of either E or I is zero. These curves provide valuable
insight into the behavior of the system. We define the nullclines using the functions I_E(E) and E_I(I), which solve the
equations dE/dt = 0 and dI/dt = 0, respectively. These functions are derived from the rate equations of the system.

12
The nullclines are plotted for a range of values of E and I, and they help visualize the locations of equilibrium points
where both derivatives are zero. The equilibria are found by solving these nullcline equations for E and I.
2.4. Equilibrium Points
Equilibrium points are the points where the system's trajectory reaches a steady state, i.e., where both dE/dt and dI/dt
are zero. These points are identified by solving the nullcline equations and locating their intersections. At these points,
the vector field is zero, and the system's behavior depends on the nature of the equilibrium (e.g., stable or unstable).
We assume that the equilibrium points are stable if the surrounding vector field tends to push the system back toward
the equilibrium, and unstable if the vector field pushes the system away.
2.5. Local Phase Portraits
Local phase portraits show the behavior of the system in the vicinity of equilibrium points. They are created by
simulating the system's trajectories for different initial conditions close to an equilibrium point and observing how the
system evolves over time. These portraits provide information on the stability of the equilibrium and the system's
dynamic behavior.
In our approach, the vector field is plotted near equilibrium points, with the vectors normalized to show the direction of
the system's evolution. The local phase portraits can be used to classify the equilibrium points as nodes, saddles, or foci,
based on the system's dynamics.
2.6. Global Phase Portraits
The global phase portrait provides a complete overview of the system's behavior over the entire state space. It is
created by simulating the system from various initial conditions and plotting the trajectories in the phase plane. This
helps in understanding the global behavior of the system, including the flow patterns, the influence of equilibrium
points, and the overall stability of the system.
The quiver plot is used to visualize the vector field, where the length and direction of the vectors indicate the rate of
change of E and I. In our implementation, the vectors are normalized and scaled to avoid overlapping or excessively long
arrows, especially near equilibrium points.
3. Results
The following results were obtained through simulation and visualization:
 Nullclines: The nullclines for E and I were computed and plotted. The curves intersect at equilibrium points,
where the system's derivatives are zero.
 Equilibrium Points: The equilibrium points were identified at the intersections of the nullclines. These points
correspond to steady states in the system where no change occurs.
 Local Phase Portraits: Local phase portraits were generated near the equilibrium points to analyze the system's
stability. The behavior of the system near equilibrium points was consistent with the expected dynamics based
on the chosen functions.
 Global Phase Portrait: A global phase portrait was constructed by simulating the system from various initial
conditions. The resulting plot showed the trajectories in the phase plane and provided a clear visual
representation of the system's behavior, including the flow patterns and the influence of equilibrium points.
The vector field was normalized and scaled to prevent overly large vectors, and its behavior near the equilibrium points
indicated that the system slows down and the vectors decrease in magnitude as expected.
4. Discussion
The results obtained make sense based on the model's structure and the chosen numerical methods. The equilibrium
points are correctly identified, and the nullclines intersect as expected. The local phase portraits show the expected
behavior near the equilibrium points, and the global phase portrait provides a clear overview of the system's dynamics.
The Euler method is adequate for this simulation, given the relatively simple nature of the system and the small time
step used. While more sophisticated methods (such as Runge-Kutta) could provide higher accuracy, the Euler method
strikes a balance between simplicity and computational efficiency for this task.
The choice of parameters, such as the time step (dt = 0.01), was made to ensure smooth trajectories without causing
excessive computational cost. The step size was small enough to accurately capture the system's dynamics while
maintaining reasonable computational time.
5. Conclusion
This study successfully simulates and visualizes the phase portrait of a two-variable dynamical system. The methods
used, including the Euler method for numerical integration and quiver plots for the vector field, were effective in
representing the system's behavior. The results obtained were consistent with expectations, and the model provides
valuable insights into the dynamics of the system.

13
NB. I have used a combination of Copilot and Chatgpt to improve my code and also work through some codes that
are relatively black box.

14