This is a matter of what model you fit and your experimental objective.
The short answer is that including a central point (a) makes it possible to detect when a more complicated model is needed and, regardless, (b) is required if you are searching for optimal responses.
Let's look at some of the possibilities.
1. Linear fit
Assuming you want to model some differential among the responses, about the simplest model that allows for (a) an average value to be estimated and (b, c) variations along both the $x$ and $y$ coordinates would take the form
$$z = \beta_0 + \beta_1 f(x) + \beta_2 g(y) + \text{error}$$
where $f$ is some specified function expected to describe variation in the $x$ coordinate (up to some multiple $\beta_1$ to be estimated) and likewise for $g$ relative to $y.$ Generally, $f$ and $g$ are "feature engineering" functions that transform the data $(x,y)$ as recorded into variables $(f(x), g(y))$ you hope are linearly related to the response $z.$ The very simplest case takes $f(x) = x$ and $g(y) = y,$ where you don't transform the data at all.
Geometrically, this model fits a plane to the data.
In this pseudo-3D image, colors indicate function values and $f$ labels the "height" or "response" axis.
It is a geometrical fact, going back to Euclid, that to determine a plane you need to make observations at a minimum of three distinct, non-collinear points in the plane. Taking the values of the features at these three points to be $(f_1,g_1),$ $(f_2,g_2),$ and $(f_3,g_3),$ an algebraic statement of this criterion is
$$\det \pmatrix {1 & f_1 & g_1 \\1 & f_2 & g_2 \\ 1 & f_3 & g_3} \ne 0.$$
This generalizes to more variables, as we will see.
2. Linear fit with interaction
This model includes one more feature formed from the product of the two features, known as an "interaction:"
$$z = \beta_0 + \beta_1 f(x) + \beta_2 g(y) + \beta_{12} f(x)g(y) + \text{error}$$
It requires one more parameter $\beta_{12}$ to be fit.
Geometrically, this model fits a ruled surface to the data.
When $f = x$ and $g = y,$ the gridlines are straight. That is because each one traces out points on the surface whether one of $x$ or $y$ is fixed. For instance, taking $y$ to be constant, any such gridline satisfies the equation
$$z = (\beta_0 + \beta_2 y) + (\beta_1 + \beta_{12}g(y)) x = \alpha_0 + \alpha_1 x,$$
which as $x$ varies describes a line.
Such a model would be indicated when the estimate of the interaction coefficient $\hat\beta_{12}$ is significant: that is, it is so much larger in size than its standard error that you would be uncomfortable not including the interaction.
A basic theorem about ruled surfaces is that every point on a smooth, non-flat ruled surface has negative Gaussian curvature. For optimal response surface design that is really bad, because at a point of negative curvature the surface curves up in one direction and down in perpendicular directions. Maxima and minima don't exist!
Incidentally, the criterion for fitting all four parameters to the data is that
$$\det \pmatrix {1 & f_1 & g_1 & f_1 g_1 \\1 & f_2 & g_2 & f_2g_2 \\ 1 & f_3 & g_3 & f_3g_3 \\ 1 & f_4 & g_4 & f_4g_4} \ne 0$$
where the subscripts $1,2,3,4$ refer to four of the points where you made measurements. (If you used more than four points, it suffices for some four of them to satisfy this criterion.)
In particular, to succeed in estimating all four parameters you need to make observations at at least four distinct points. When the points are placed at vertices of a rectangle parallel to the coordinate axes, as in the question, this determinant will be nonzero, assuring a unique solution.
3. Fits with positive curvature
At any smooth, non-flat local extremum, the Gaussian curvature is positive. In order to create that possibility in the model, at least one more term must be included. Ordinarily we would include terms for both $f(x)^2$ and $g(y)^2,$ but if we wish to be truly parsimonious we might suppose their coefficients are the same ("isotropy"), so that only one more parameter is needed. Now the model is
$$z = \beta_0 + \beta_1 f(x) + \beta_2 g(y) + \beta_{12} f(x)g(y) + \beta_3(f(x)^2 + g(y)^2) + \text{error}$$
For the first time, the surface may exhibit an extremum:
However, in order to identify this surface uniquely, observations must be made on at least five distinct locations. The criterion for successfully fitting all five parameters is now
$$\det \pmatrix {1 & f_1 & g_1 & f_1 g_1 & f_1^2 + g_1^2\\1 & f_2 & g_2 & f_2g_2 & f_2^2 + g_2^2\\ 1 & f_3 & g_3 & f_3g_3 & f_3^2 + g_3^2\\ 1 & f_4 & g_4 & f_4g_4& f_4^2 + g_4^2\\ 1 & f_5 & g_5 & f_5g_5& f_5^2 + g_5^2} \ne 0$$
for some five points of the experiment. This will be the case for the center-point design shown in the question.
Even if this new coefficient estimate $\hat\beta_3$ appears non-significant, you might want to include it anyway. After all, without it there is no possibility of finding a local optimum within the interior (convex hull) of your design points. Moreover, you shouldn't be testing $\beta_3$ by itself for significance. It's really part of the two sets of factors (the interaction and quadratic term) the model departures from the initial linear fit. It's possible for each of $\hat\beta_{12}$ and $\hat\beta_3$ to be individually not significant but for the two of them together to be significant.
